diff --git a/pytential/muller.py b/pytential/muller.py
index 22fe5458e895236e3ec0fcc6ab7634290e0a84ed..5fc6a82378af0e8c86ac5858187014399cf8c42d 100644
--- a/pytential/muller.py
+++ b/pytential/muller.py
@@ -25,7 +25,7 @@ THE SOFTWARE.
 import numpy as np
 
 
-def muller_deflate(f, n, maxiter=100, eps=1e-14):
+def muller_deflate(f, n, maxiter=100, eps=1e-14, z_start=None):
     """
     :arg n: number of zeros sought
     :returns: (roots, niter) - roots of the given function; number of
@@ -52,7 +52,7 @@ def muller_deflate(f, n, maxiter=100, eps=1e-14):
         niter.append(niter0)
 
         while (np.isnan(roots[i]) or niter[i] == maxiter) and miter < 50:
-            roots0, niter0 = muller(f_deflated, maxiter, eps)
+            roots0, niter0 = muller(f_deflated, maxiter, eps, z_start=z_start)
             roots[i] = roots0
             niter[i] = niter0
             miter = miter+1
diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index 9c9b4ea399f56895c10b80003752c56f6fe95e98..0d058c820944180ed27bea27579051cb71fc9b5b 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -25,7 +25,6 @@ THE SOFTWARE.
 
 import six
 
-import loopy as lp
 import numpy as np
 from pytools import memoize_method
 from meshmode.discretization import Discretization
@@ -45,45 +44,6 @@ __doc__ = """
 """
 
 
-# {{{ jump term interface helper
-
-class _JumpTermArgumentProvider(object):
-    def __init__(self, discr, density, ds_direction, side=None):
-        self.discr = discr
-        self.density = density
-        self.ds_direction = ds_direction
-        self.side = side
-
-    @property
-    def normal(self):
-        return self.discr.curve.normals.reshape(2, -1).T
-
-    @property
-    def tangent(self):
-        return self.discr.curve.tangents.reshape(2, -1).T
-
-    @property
-    def src_derivative_dir(self):
-        return self.ds_direction
-
-    @property
-    def mean_curvature(self):
-        return self.discr.curve.curvature.reshape(-1)
-
-    @property
-    def density_0(self):
-        return self.density.reshape(-1)
-
-    @property
-    @memoize_method
-    def density_0_prime(self):
-        diff_mat = self.discr.curve.expansion.get_differentiation_matrix()
-        return (2 * np.dot(diff_mat, self.density.T).T.reshape(-1)
-                / self.discr.curve.speed.reshape(-1))
-
-# }}}
-
-
 class LayerPotentialSource(PotentialSource):
     pass
 
@@ -128,13 +88,14 @@ class QBXLayerPotentialSource(LayerPotentialSource):
     .. attribute :: qbx_order
     .. attribute :: fmm_order
     .. attribute :: cl_context
-    .. automethod :: centers
-    .. automethod :: panel_sizes
     .. automethod :: weights_and_area_elements
     .. automethod :: with_refinement
 
     See :ref:`qbxguts` for some information on the inner workings of this.
     """
+
+    # {{{ constructor / copy
+
     def __init__(self, density_discr, fine_order,
             qbx_order=None,
             fmm_order=None,
@@ -146,6 +107,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
             # FIXME default debug=False once everything works
             debug=True,
             refined_for_global_qbx=False,
+            expansion_disks_in_tree_have_extent=False,
             performance_data_file=None):
         """
         :arg fine_order: The total degree to which the (upsampled)
@@ -184,6 +146,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
 
         self.debug = debug
         self.refined_for_global_qbx = refined_for_global_qbx
+        self.expansion_disks_in_tree_have_extent = \
+                expansion_disks_in_tree_have_extent
         self.performance_data_file = performance_data_file
 
     def copy(
@@ -216,12 +180,17 @@ class QBXLayerPotentialSource(LayerPotentialSource):
                 refined_for_global_qbx=(
                     refined_for_global_qbx if refined_for_global_qbx is not None
                     else self.refined_for_global_qbx),
+                expansion_disks_in_tree_have_extent=(
+                    self.expansion_disks_in_tree_have_extent),
                 performance_data_file=self.performance_data_file)
 
+    # }}}
+
     @property
     def base_fine_density_discr(self):
-        """The fine density discretization before upsampling is applied.
+        """The refined, interpolation-focused density discretization (no oversampling).
         """
+        # FIXME: Maybe rename refined_interp_density_discr
         return (self._base_resampler.to_discr
                 if self._base_resampler is not None
                 else self.density_discr)
@@ -229,6 +198,9 @@ class QBXLayerPotentialSource(LayerPotentialSource):
     @property
     @memoize_method
     def fine_density_discr(self):
+        """The refined, quadrature-focused density discretization (with upsampling).
+        """
+        # FIXME: Maybe rename refined_quad_density_discr
         from meshmode.discretization.poly_element import (
                 QuadratureSimplexGroupFactory)
 
@@ -251,22 +223,6 @@ class QBXLayerPotentialSource(LayerPotentialSource):
 
         return conn
 
-    def el_view(self, discr, group_nr, global_array):
-        """Return a view of *global_array* of shape
-        ``(..., discr.groups[group_nr].nelements)``
-        where *global_array* is of shape ``(..., nelements)``,
-        where *nelements* is the global (per-discretization) node count.
-        """
-
-        group = discr.groups[group_nr]
-        el_nr_base = sum(group.nelements for group in discr.groups[:group_nr])
-
-        return global_array[
-            ..., el_nr_base:el_nr_base + group.nelements] \
-            .reshape(
-                global_array.shape[:-1]
-                + (group.nelements,))
-
     @property
     @memoize_method
     def refiner_code_container(self):
@@ -305,7 +261,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
     @memoize_method
     def h_max(self):
         with cl.CommandQueue(self.cl_context) as queue:
-            panel_sizes = self.panel_sizes("npanels").with_queue(queue)
+            panel_sizes = self._panel_sizes("npanels").with_queue(queue)
             return np.asscalar(cl.array.max(panel_sizes).get())
 
     @property
@@ -328,121 +284,72 @@ class QBXLayerPotentialSource(LayerPotentialSource):
     def complex_dtype(self):
         return self.density_discr.complex_dtype
 
-    def _centers_of_mass_for_discr(self, discr):
-        knl = lp.make_kernel(
-            """{[dim,k,i]:
-                0<=dim<ndims and
-                0<=k<nelements and
-                0<=i<nunit_nodes}""",
-            """
-                panels[dim, k] = sum(i, nodes[dim, k, i])/nunit_nodes
-                """,
-            default_offset=lp.auto, name="find_panel_centers_of_mass")
-
-        knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
-        knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
-        knl = lp.tag_inames(knl, dict(dim="ilp"))
-
-        with cl.CommandQueue(self.cl_context) as queue:
-            mesh = discr.mesh
-            panels = cl.array.empty(queue, (mesh.ambient_dim, mesh.nelements),
-                                    dtype=discr.real_dtype)
-            for group_nr, group in enumerate(discr.groups):
-                _, (result,) = knl(queue,
-                    nelements=group.nelements,
-                    nunit_nodes=group.nunit_nodes,
-                    nodes=group.view(discr.nodes()),
-                    panels=self.el_view(discr, group_nr, panels))
-            panels.finish()
-            panels = panels.with_queue(None)
-            return tuple(panels[d, :] for d in range(mesh.ambient_dim))
+    # {{{ internal API
 
     @memoize_method
-    def panel_centers_of_mass(self):
-        return self._centers_of_mass_for_discr(self.density_discr)
+    def _panel_centers_of_mass(self):
+        import pytential.qbx.utils as utils
+        return utils.element_centers_of_mass(self.density_discr)
 
     @memoize_method
-    def fine_panel_centers_of_mass(self):
-        return self._centers_of_mass_for_discr(self.base_fine_density_discr)
-
-    def _panel_sizes_for_discr(self, discr, last_dim_length):
-        assert last_dim_length in ("nsources", "ncenters", "npanels")
-        # To get the panel size this does the equivalent of (∫ 1 ds)**(1/dim).
-        # FIXME: Kernel optimizations
+    def _fine_panel_centers_of_mass(self):
+        import pytential.qbx.utils as utils
+        return utils.element_centers_of_mass(self.base_fine_density_discr)
 
-        if last_dim_length == "nsources" or last_dim_length == "ncenters":
-            knl = lp.make_kernel(
-                "{[i,j,k]: 0<=i<nelements and 0<=j,k<nunit_nodes}",
-                "panel_sizes[i,j] = sum(k, ds[i,k])**(1/dim)",
-                name="compute_size")
-
-            def panel_size_view(discr, group_nr):
-                return discr.groups[group_nr].view
-
-        elif last_dim_length == "npanels":
-            knl = lp.make_kernel(
-                "{[i,j]: 0<=i<nelements and 0<=j<nunit_nodes}",
-                "panel_sizes[i] = sum(j, ds[i,j])**(1/dim)",
-                name="compute_size")
-            from functools import partial
+    @memoize_method
+    def _expansion_radii(self, last_dim_length):
+        if last_dim_length == "npanels":
+            # FIXME: Make this an error
 
-            def panel_size_view(discr, group_nr):
-                return partial(self.el_view, discr, group_nr)
+            from warnings import warn
+            warn("Passing 'npanels' as last_dim_length to _expansion_radii is "
+                    "deprecated. Expansion radii should be allowed to vary "
+                    "within a panel.", stacklevel=3)
 
-        else:
-            raise ValueError("unknown dim length specified")
+        with cl.CommandQueue(self.cl_context) as queue:
+                return (self._panel_sizes(last_dim_length).with_queue(queue) * 0.5
+                        ).with_queue(None)
 
-        knl = lp.fix_parameters(knl, dim=self.dim)
+    # _expansion_radii should not be needed for the fine discretization
 
+    @memoize_method
+    def _close_target_tunnel_radius(self, last_dim_length):
         with cl.CommandQueue(self.cl_context) as queue:
-            from pytential import bind, sym
-            ds = bind(
-                    discr,
-                    sym.area_element(ambient_dim=discr.ambient_dim, dim=discr.dim)
-                    * sym.QWeight()
-                    )(queue)
-            panel_sizes = cl.array.empty(
-                queue, discr.nnodes
-                if last_dim_length in ("nsources", "ncenters")
-                else discr.mesh.nelements, discr.real_dtype)
-            for group_nr, group in enumerate(discr.groups):
-                _, (result,) = knl(queue,
-                    nelements=group.nelements,
-                    nunit_nodes=group.nunit_nodes,
-                    ds=group.view(ds),
-                    panel_sizes=panel_size_view(
-                        discr, group_nr)(panel_sizes))
-            panel_sizes.finish()
-            if last_dim_length == "ncenters":
-                from pytential.qbx.utils import get_interleaver_kernel
-                knl = get_interleaver_kernel(discr.real_dtype)
-                _, (panel_sizes,) = knl(queue, dstlen=2*discr.nnodes,
-                                        src1=panel_sizes, src2=panel_sizes)
-            return panel_sizes.with_queue(None)
+                return (self._panel_sizes(last_dim_length).with_queue(queue) * 0.5
+                        ).with_queue(None)
 
     @memoize_method
-    def panel_sizes(self, last_dim_length="npanels"):
-        return self._panel_sizes_for_discr(self.density_discr, last_dim_length)
+    def _panel_sizes(self, last_dim_length="npanels"):
+        import pytential.qbx.utils as utils
+        return utils.panel_sizes(self.density_discr, last_dim_length)
 
     @memoize_method
-    def fine_panel_sizes(self, last_dim_length="npanels"):
+    def _fine_panel_sizes(self, last_dim_length="npanels"):
         if last_dim_length != "npanels":
             raise NotImplementedError()
-        return self._panel_sizes_for_discr(
-                self.base_fine_density_discr, last_dim_length)
+
+        import pytential.qbx.utils as utils
+        return utils.panel_sizes(self.base_fine_density_discr, last_dim_length)
 
     @memoize_method
-    def centers(self, sign):
-        adim = self.density_discr.ambient_dim
-        dim = self.density_discr.dim
+    def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides):
+        """
+        :arg target_discrs_and_qbx_sides:
+            a tuple of *(discr, qbx_forced_limit)*
+            tuples, where *discr* is a
+            :class:`meshmode.discretization.Discretization`
+            or
+            :class:`pytential.target.TargetBase`
+            instance
+        """
+        from pytential.qbx.geometry import QBXFMMGeometryData
 
-        from pytential import sym, bind
-        with cl.CommandQueue(self.cl_context) as queue:
-            nodes = bind(self.density_discr, sym.nodes(adim))(queue)
-            normals = bind(self.density_discr, sym.normal(adim, dim=dim))(queue)
-            panel_sizes = self.panel_sizes("nsources").with_queue(queue)
-            return (nodes + normals * sign * panel_sizes / 2).as_vector(np.object)
+        return QBXFMMGeometryData(self.qbx_fmm_code_getter,
+                self, target_discrs_and_qbx_sides,
+                target_stick_out_factor=self.target_stick_out_factor,
+                debug=self.debug)
+
+    # }}}
 
     @memoize_method
     def weights_and_area_elements(self):
@@ -463,6 +370,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
 
             return (area_element.with_queue(queue)*qweight).with_queue(None)
 
+    # {{{ helpers for symbolic operator processing
+
     def preprocess_optemplate(self, name, discretizations, expr):
         """
         :arg name: The symbolic name for *self*, which the preprocessor
@@ -480,6 +389,10 @@ class QBXLayerPotentialSource(LayerPotentialSource):
 
         return result
 
+    # }}}
+
+    # {{{ internal functionality for execution
+
     def exec_compute_potential_insn(self, queue, insn, bound_expr, evaluate):
         from pytools.obj_array import with_object_array_or_scalar
         from functools import partial
@@ -509,24 +422,6 @@ class QBXLayerPotentialSource(LayerPotentialSource):
         return QBXFMMGeometryCodeGetter(self.cl_context, self.ambient_dim,
                 debug=self.debug)
 
-    @memoize_method
-    def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides):
-        """
-        :arg target_discrs_and_qbx_sides:
-            a tuple of *(discr, qbx_forced_limit)*
-            tuples, where *discr* is a
-            :class:`meshmode.discretization.Discretization`
-            or
-            :class:`pytential.target.TargetBase`
-            instance
-        """
-        from pytential.qbx.geometry import QBXFMMGeometryData
-
-        return QBXFMMGeometryData(self.qbx_fmm_code_getter,
-                self, target_discrs_and_qbx_sides,
-                target_stick_out_factor=self.target_stick_out_factor,
-                debug=self.debug)
-
     # {{{ fmm-based execution
 
     @memoize_method
@@ -650,7 +545,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
 
         # }}}
 
-        if len(geo_data.global_qbx_centers()) != geo_data.center_info().ncenters:
+        if len(geo_data.global_qbx_centers()) != geo_data.ncenters:
             raise NotImplementedError("geometry has centers requiring local QBX")
 
         from pytential.qbx.geometry import target_state
@@ -768,6 +663,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
         strengths = (evaluate(insn.density).with_queue(queue)
                 * self.weights_and_area_elements())
 
+        import pytential.qbx.utils as utils
+
         # FIXME: Do this all at once
         result = []
         for o in insn.outputs:
@@ -783,7 +680,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
                 evt, output_for_each_kernel = lpot_applier(
                         queue, target_discr.nodes(),
                         self.fine_density_discr.nodes(),
-                        self.centers(o.qbx_forced_limit),
+                        utils.get_centers_on_side(self, o.qbx_forced_limit),
                         [strengths], **kernel_args)
                 result.append((o.name, output_for_each_kernel[o.kernel_index]))
             else:
@@ -807,12 +704,11 @@ class QBXLayerPotentialSource(LayerPotentialSource):
                             (target_discr, qbx_forced_limit),
                         ))
 
-                # center_info is independent of targets
-                center_info = geo_data.center_info()
+                # center-related info is independent of targets
 
                 # First ncenters targets are the centers
                 tgt_to_qbx_center = (
-                        geo_data.user_target_to_center()[center_info.ncenters:]
+                        geo_data.user_target_to_center()[geo_data.ncenters:]
                         .copy(queue=queue))
 
                 qbx_tgt_numberer = self.get_qbx_target_numberer(
@@ -845,7 +741,7 @@ class QBXLayerPotentialSource(LayerPotentialSource):
                             queue,
                             targets=target_discr.nodes(),
                             sources=self.fine_density_discr.nodes(),
-                            centers=center_info.centers,
+                            centers=geo_data.centers(),
                             strengths=[strengths],
                             qbx_tgt_numbers=qbx_tgt_numbers,
                             qbx_center_numbers=qbx_center_numbers,
@@ -857,6 +753,8 @@ class QBXLayerPotentialSource(LayerPotentialSource):
 
     # }}}
 
+    # }}}
+
 # }}}
 
 
diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py
index e3508b7089d2875c12a77079e763201300b5a30a..7fd7f49d4e957d216d0fd83c680740fb11e40521 100644
--- a/pytential/qbx/fmm.py
+++ b/pytential/qbx/fmm.py
@@ -153,7 +153,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
         return cl.array.zeros(
                     self.queue,
-                    (self.geo_data.center_info().ncenters,
+                    (self.geo_data.ncenters,
                         len(qbx_l_expn)),
                     dtype=self.dtype)
 
@@ -212,7 +212,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
                 self.queue,
                 global_qbx_centers=geo_data.global_qbx_centers(),
                 qbx_center_to_target_box=geo_data.qbx_center_to_target_box(),
-                qbx_centers=geo_data.center_info().centers,
+                qbx_centers=geo_data.centers(),
 
                 source_box_starts=starts,
                 source_box_lists=lists,
@@ -230,7 +230,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
         qbx_expansions = self.qbx_local_expansion_zeros()
 
         geo_data = self.geo_data
-        if geo_data.center_info().ncenters == 0:
+        if geo_data.ncenters == 0:
             return qbx_expansions
 
         traversal = geo_data.traversal()
@@ -249,7 +249,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
                     qbx_center_to_target_box=geo_data.qbx_center_to_target_box(),
 
                     centers=self.tree.box_centers,
-                    qbx_centers=geo_data.center_info().centers,
+                    qbx_centers=geo_data.centers(),
 
                     src_expansions=source_mpoles_view,
                     src_base_ibox=source_level_start_ibox,
@@ -273,7 +273,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
         qbx_expansions = self.qbx_local_expansion_zeros()
 
         geo_data = self.geo_data
-        if geo_data.center_info().ncenters == 0:
+        if geo_data.ncenters == 0:
             return qbx_expansions
         trav = geo_data.traversal()
 
@@ -294,7 +294,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
                     target_base_ibox=target_level_start_ibox,
 
                     centers=self.tree.box_centers,
-                    qbx_centers=geo_data.center_info().centers,
+                    qbx_centers=geo_data.centers(),
 
                     expansions=target_locals_view,
                     qbx_expansions=qbx_expansions,
@@ -322,7 +322,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
         qbxl2p = self.code.qbxl2p(self.qbx_order)
 
         evt, pot_res = qbxl2p(self.queue,
-                qbx_centers=geo_data.center_info().centers,
+                qbx_centers=geo_data.centers(),
                 global_qbx_centers=geo_data.global_qbx_centers(),
 
                 center_to_targets_starts=ctt.starts,
@@ -652,7 +652,7 @@ def write_performance_model(outf, geo_data):
     qbx_center_to_target_box = geo_data.qbx_center_to_target_box()
     center_to_targets_starts = geo_data.center_to_tree_targets().starts
 
-    ncenters = geo_data.center_info().ncenters
+    ncenters = geo_data.ncenters
 
     outf.write("ncenters = {cost}\n"
             .format(cost=ncenters))
diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py
index ce93e851d61741af3951096d3f8e131c76b221fc..0eee50d65f1f8ea61cd8f4f960f14482abd6b464 100644
--- a/pytential/qbx/geometry.py
+++ b/pytential/qbx/geometry.py
@@ -118,72 +118,6 @@ class QBXFMMGeometryCodeGetter(object):
         self.ambient_dim = ambient_dim
         self.debug = debug
 
-    @property
-    @memoize_method
-    def pick_expansion_centers(self):
-        knl = lp.make_kernel(
-            """{[dim,k,i]:
-                0<=dim<ndims and
-                0<=k<nelements and
-                0<=i<nout_nodes}""",
-            """
-                centers[dim, k, i] = all_centers[dim, k, kept_center_indices[i]]
-                radii[k, i] = all_radii[k, kept_center_indices[i]]
-                """,
-            [
-                lp.GlobalArg("all_centers", None,
-                    shape="ndims,nelements,nunit_nodes"),
-                lp.GlobalArg("all_radii", None, shape="nelements,nunit_nodes"),
-                lp.ValueArg("nunit_nodes", np.int32),
-                "..."
-                ],
-            default_offset=lp.auto,
-            name="center_pick")
-
-        knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
-        knl = lp.tag_array_axes(knl, "centers,all_centers", "sep, C, C")
-
-        knl = lp.split_iname(knl, "i", 16, inner_tag="l.0")
-        return lp.tag_inames(knl, dict(k="g.0", dim="ilp"))
-
-    @property
-    @memoize_method
-    def find_element_centers(self):
-        knl = lp.make_kernel(
-            """{[dim,k,i]:
-                0<=dim<ndims and
-                0<=k<nelements and
-                0<=i<nunit_nodes}""",
-            """
-                el_centers[dim, k] = sum(i, nodes[dim, k, i])/nunit_nodes
-                """,
-            default_offset=lp.auto, name="find_element_centers")
-
-        knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
-        knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
-        return lp.tag_inames(knl, dict(dim="ilp"))
-
-    @property
-    @memoize_method
-    def find_element_radii(self):
-        knl = lp.make_kernel(
-            """{[dim,k,i]:
-                0<=dim<ndims and
-                0<=k<nelements and
-                0<=i<nunit_nodes}""",
-            """
-                el_radii[k] = max(dim, max(i, \
-                    fabs(nodes[dim, k, i] - el_centers[dim, k])))
-                """,
-            default_offset=lp.auto, name="find_element_radii")
-
-        knl = lp.fix_parameters(knl, ndims=self.ambient_dim)
-
-        knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
-        return lp.tag_inames(knl, dict(dim="unr"))
-
     @memoize_method
     def copy_targets_kernel(self):
         knl = lp.make_kernel(
@@ -263,104 +197,6 @@ class QBXFMMGeometryCodeGetter(object):
         from boxtree.area_query import LeavesToBallsLookupBuilder
         return LeavesToBallsLookupBuilder(self.cl_context)
 
-    # {{{ check if a center may be used with global QBX
-
-    @memoize_method
-    def qbx_center_for_global_tester(self,
-            coord_dtype, box_id_dtype, particle_id_dtype):
-        from pyopencl.elementwise import ElementwiseTemplate
-        return ElementwiseTemplate(
-            arguments="""//CL:mako//
-                /* input */
-                %for iaxis in range(ambient_dim):
-                    coord_t *center_${iaxis},  /* [ncenters] */
-                %endfor
-                coord_t *radii,  /* [ncenters] */
-                box_id_t *qbx_center_to_target_box, /* [ncenters] */
-
-                box_id_t *neighbor_source_boxes_starts,
-                box_id_t *neighbor_source_boxes_lists,
-                particle_id_t *box_point_source_starts,
-                particle_id_t *box_point_source_counts_cumul,
-                %for iaxis in range(ambient_dim):
-                    coord_t *point_sources_${iaxis},
-                %endfor
-
-                /* output */
-                char *center_may_use_global_qbx,
-                """,
-            operation=r"""//CL:mako//
-                particle_id_t icenter = i;
-
-                %for iaxis in range(ambient_dim):
-                    coord_t my_center_${iaxis} = center_${iaxis}[icenter];
-                %endfor
-                coord_t radius = radii[icenter];
-
-                // {{{ see if there are sources close enough to require local QBX
-
-                bool found_too_close = false;
-
-                coord_t radius_squared = radius * radius;
-
-                box_id_t itgt_box = qbx_center_to_target_box[icenter];
-
-                box_id_t nb_src_start = neighbor_source_boxes_starts[itgt_box];
-                box_id_t nb_src_stop = neighbor_source_boxes_starts[itgt_box+1];
-
-                for (
-                        box_id_t ibox_list = nb_src_start;
-                        ibox_list < nb_src_stop && !found_too_close;
-                        ++ibox_list)
-                {
-                    box_id_t neighbor_box_id =
-                        neighbor_source_boxes_lists[ibox_list];
-
-                    box_id_t bps_start = box_point_source_starts[neighbor_box_id];
-                    box_id_t bps_stop =
-                        bps_start + box_point_source_counts_cumul[neighbor_box_id];
-                    for (
-                            box_id_t ipsrc = bps_start;
-                            ipsrc < bps_stop;
-                            ++ipsrc)
-                    {
-                        %for iaxis in range(ambient_dim):
-                            coord_t psource_${iaxis} = point_sources_${iaxis}[ipsrc];
-                        %endfor
-
-                        coord_t dist_squared = 0
-                        %for iaxis in range(ambient_dim):
-                            + (my_center_${iaxis} - psource_${iaxis})
-                              * (my_center_${iaxis} - psource_${iaxis})
-                        %endfor
-                            ;
-
-                        const coord_t too_close_slack = (1+1e-2)*(1+1e-2);
-                        found_too_close = found_too_close
-                            || (dist_squared*too_close_slack  < radius_squared);
-
-                        if (found_too_close)
-                            break;
-                    }
-                }
-
-                // }}}
-
-                center_may_use_global_qbx[icenter] = !found_too_close;
-                """,
-            name="qbx_test_center_for_global").build(
-                    self.cl_context,
-                    type_aliases=(
-                        ("box_id_t", box_id_dtype),
-                        ("particle_id_t", particle_id_dtype),
-                        ("coord_t", coord_dtype),
-                        ),
-                    var_values=(
-                        ("ambient_dim", self.ambient_dim),
-                        ))
-
-    # }}}
-
     @property
     @memoize_method
     def key_value_sort(self):
@@ -379,7 +215,10 @@ class QBXFMMGeometryCodeGetter(object):
                     VectorArg(particle_id_dtype, "filtered_target_id"),
                     VectorArg(particle_id_dtype, "count"),
                     ],
+
+                # "Does this target have a QBX center?"
                 input_expr="(target_to_center[i] >= 0) ? 1 : 0",
+
                 scan_expr="a+b", neutral="0",
                 output_statement="""
                     if (prev_item != item)
@@ -419,32 +258,6 @@ class TargetInfo(DeviceDataRecord):
     """
 
 
-class CenterInfo(DeviceDataRecord):
-    """Information on location of QBX centers.
-    Returned from :meth:`QBXFMMGeometryData.center_info`.
-
-    .. attribute:: centers
-
-        Shape: ``[dim][ncenters]``
-
-    .. attribute:: sides
-
-        Shape: ``[ncenters]``
-
-        -1 for inside, +1 for outside (relative to normal)
-
-    .. attribute:: radii
-
-        Shape: ``[ncenters]``
-
-    .. attribute:: ncenters
-    """
-
-    @property
-    def ncenters(self):
-        return len(self.radii)
-
-
 class CenterToTargetList(DeviceDataRecord):
     """A lookup table of targets covered by each QBX disk. Indexed by global
     number of QBX center, ``lists[start[i]:start[i+1]]`` indicates numbers
@@ -516,19 +329,22 @@ class QBXFMMGeometryData(object):
 
     .. attribute:: coord_dtype
 
+    .. rubric :: Expansion centers
+
+    .. attribute:: ncenters
+    .. automethod:: centers()
+    .. automethod:: radii()
+
     .. rubric :: Methods
 
-    .. automethod:: center_info()
     .. automethod:: target_info()
     .. automethod:: tree()
     .. automethod:: traversal()
-    .. automethod:: leaf_to_center_lookup
     .. automethod:: qbx_center_to_target_box()
     .. automethod:: global_qbx_flags()
     .. automethod:: global_qbx_centers()
     .. automethod:: user_target_to_center()
     .. automethod:: center_to_tree_targets()
-    .. automethod:: global_qbx_centers_box_target_lists()
     .. automethod:: non_qbx_box_target_lists()
     .. automethod:: plot()
     """
@@ -561,35 +377,36 @@ class QBXFMMGeometryData(object):
     def coord_dtype(self):
         return self.lpot_source.fine_density_discr.nodes().dtype
 
-    # {{{ center info
+    # {{{ centers/radii
 
-    @memoize_method
-    def kept_center_indices(self, el_group):
-        """Return indices of unit nodes (of the target discretization)
-        that will give rise to QBX centers.
-        """
-
-        # FIXME: Be more careful about which nodes to keep
-        return np.arange(0, el_group.nunit_nodes)
+    @property
+    def ncenters(self):
+        return len(self.centers()[0])
 
     @memoize_method
-    def center_info(self):
-        """ Return a :class:`CenterInfo`. |cached|
-        """
+    def centers(self):
+        """ Return an object array of (interleaved) center coordinates.
 
-        lpot_source = self.lpot_source
+        ``coord_t [ambient_dim][ncenters]``
+        """
 
         with cl.CommandQueue(self.cl_context) as queue:
             from pytential.qbx.utils import get_interleaved_centers
-            centers = get_interleaved_centers(queue, lpot_source)
-            sides = cl.array.arange(queue, len(centers[0]), dtype=np.int32)
-            sides = 2 * (sides & 1) - 1
-            radii = lpot_source.panel_sizes("ncenters").with_queue(queue) / 2
+            from pytools.obj_array import make_obj_array
+            return make_obj_array([
+                ccomp.with_queue(None)
+                for ccomp in get_interleaved_centers(queue, self.lpot_source)])
 
-        return CenterInfo(
-                sides=sides,
-                radii=radii,
-                centers=centers).with_queue(None)
+    @memoize_method
+    def expansion_radii(self):
+        """Return an  array of radii associated with the (interleaved)
+        expansion centers.
+
+        ``coord_t [ncenters]``
+        """
+        with cl.CommandQueue(self.cl_context) as queue:
+            from pytential.qbx.utils import get_interleaved_radii
+            return get_interleaved_radii(queue, self.lpot_source)
 
     # }}}
 
@@ -603,9 +420,7 @@ class QBXFMMGeometryData(object):
         lpot_src = self.lpot_source
 
         with cl.CommandQueue(self.cl_context) as queue:
-            center_info = self.center_info()
-
-            ntargets = center_info.ncenters
+            ntargets = self.ncenters
             target_discr_starts = []
 
             for target_discr, qbx_side in self.target_discrs_and_qbx_sides:
@@ -619,8 +434,8 @@ class QBXFMMGeometryData(object):
                     self.coord_dtype)
             code_getter.copy_targets_kernel()(
                     queue,
-                    targets=targets[:, :center_info.ncenters],
-                    points=center_info.centers)
+                    targets=targets[:, :self.ncenters],
+                    points=self.centers())
 
             for start, (target_discr, _) in zip(
                     target_discr_starts, self.target_discrs_and_qbx_sides):
@@ -634,23 +449,6 @@ class QBXFMMGeometryData(object):
                     target_discr_starts=target_discr_starts,
                     ntargets=ntargets).with_queue(None)
 
-    def target_radii(self):
-        """Shape: ``[ntargets]``
-
-        A list of target radii for FMM tree construction. In this list, the QBX
-        centers have the radii determined by :meth:`center_info`, and all other
-        targets have radius zero.
-        """
-
-        tgt_info = self.target_info()
-        center_info = self.center_info()
-
-        with cl.CommandQueue(self.cl_context) as queue:
-            target_radii = cl.array.zeros(queue, tgt_info.ntargets, self.coord_dtype)
-            target_radii[:center_info.ncenters] = center_info.radii
-
-            return target_radii.with_queue(None)
-
     def target_side_preferences(self):
         """Return one big array combining all the data from
         the *side* part of :attr:`TargetInfo.target_discrs_and_qbx_sides`.
@@ -658,12 +456,11 @@ class QBXFMMGeometryData(object):
         Shape: ``[ntargets]``, dtype: int8"""
 
         tgt_info = self.target_info()
-        center_info = self.center_info()
 
         with cl.CommandQueue(self.cl_context) as queue:
             target_side_preferences = cl.array.empty(
                     queue, tgt_info.ntargets, np.int8)
-            target_side_preferences[:center_info.ncenters] = 0
+            target_side_preferences[:self.ncenters] = 0
 
             for tdstart, (target_discr, qbx_side) in \
                     zip(tgt_info.target_discr_starts,
@@ -679,7 +476,7 @@ class QBXFMMGeometryData(object):
 
     @memoize_method
     def tree(self):
-        """Build and return a :class:`boxtree.tree.TreeWithLinkedPointSources`
+        """Build and return a :class:`boxtree.tree.Tree`
         for this source with these targets.
 
         |cached|
@@ -693,6 +490,12 @@ class QBXFMMGeometryData(object):
             nsources = lpot_src.fine_density_discr.nnodes
             nparticles = nsources + target_info.ntargets
 
+            target_radii = None
+            if self.lpot_source.expansion_disks_in_tree_have_extent:
+                target_radii = cl.array.zeros(queue, target_info.ntargets,
+                        self.coord_dtype)
+                target_radii[:self.ncenters] = self.expansion_radii()
+
             refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32)
             refine_weights[:nsources] = 1
             refine_weights.finish()
@@ -707,13 +510,16 @@ class QBXFMMGeometryData(object):
             # to avoid having too many disks overlap more than one box.
             #
             # FIXME: Should investigate this further.
+
             tree, _ = code_getter.build_tree(queue,
                     particles=lpot_src.fine_density_discr.nodes(),
                     targets=target_info.targets,
+                    target_radii=target_radii,
                     max_leaf_refine_weight=256,
                     refine_weights=refine_weights,
                     debug=self.debug,
-                    kind="adaptive-level-restricted")
+                    stick_out_factor=0,
+                    kind="adaptive")
 
             return tree
 
@@ -731,31 +537,21 @@ class QBXFMMGeometryData(object):
             trav, _ = self.code_getter.build_traversal(queue, self.tree(),
                     debug=self.debug)
 
-            return trav
+            if self.lpot_source.expansion_disks_in_tree_have_extent:
+                trav = trav.merge_close_lists(queue)
 
-    def leaf_to_center_lookup(self):
-        """Return a :class:`boxtree.area_query.LeavesToBallsLookup` to look up
-        which which QBX disks overlap each leaf box.
-        """
-
-        center_info = self.center_info()
-
-        with cl.CommandQueue(self.cl_context) as queue:
-            lbl, _ = self.code_getter.build_leaf_to_ball_lookup(queue,
-                    self.tree(), center_info.centers, center_info.radii)
-            return lbl
+            return trav
 
     @memoize_method
     def qbx_center_to_target_box(self):
-        """Return a lookup table of length :attr:`CenterInfo.ncenters`
-        (see :meth:`center_info`) indicating the target box in which each
+        """Return a lookup table of length :attr:`ncenters`
+        indicating the target box in which each
         QBX disk is located.
 
         |cached|
         """
         tree = self.tree()
         trav = self.traversal()
-        center_info = self.center_info()
 
         qbx_center_to_target_box_lookup = \
                 self.code_getter.qbx_center_to_target_box_lookup(
@@ -788,18 +584,17 @@ class QBXFMMGeometryData(object):
                     box_target_starts=tree.box_target_starts,
                     box_target_counts_nonchild=tree.box_target_counts_nonchild,
                     user_target_from_tree_target=user_target_from_tree_target,
-                    ncenters=center_info.ncenters)
+                    ncenters=self.ncenters)
 
             return qbx_center_to_target_box.with_queue(None)
 
     @memoize_method
     def global_qbx_flags(self):
         """Return an array of :class:`numpy.int8` of length
-        :attr:`CenterInfo.ncenters` (see :meth:`center_info`) indicating
-        whether each center can use gloal QBX, i.e. whether a single expansion
-        can mediate interactions from *all* sources to all targets for which it
-        is valid. If global QBX can be used, the center's entry will be 1,
-        otherwise it will be 0.
+        :attr:`ncenters` indicating whether each center can use gloal QBX, i.e.
+        whether a single expansion can mediate interactions from *all* sources
+        to all targets for which it is valid. If global QBX can be used, the
+        center's entry will be 1, otherwise it will be 0.
 
         (If not, local QBX is needed, and the center may only be
         able to mediate some of the interactions to a given target.)
@@ -807,10 +602,8 @@ class QBXFMMGeometryData(object):
         |cached|
         """
 
-        center_info = self.center_info()
-
         with cl.CommandQueue(self.cl_context) as queue:
-            result = cl.array.empty(queue, center_info.ncenters, np.int8)
+            result = cl.array.empty(queue, self.ncenters, np.int8)
             result.fill(1)
 
         return result.with_queue(None)
@@ -829,7 +622,7 @@ class QBXFMMGeometryData(object):
 
             logger.info("find global qbx centers: start")
             result, count, _ = copy_if(
-                    cl.array.arange(queue, self.center_info().ncenters,
+                    cl.array.arange(queue, self.ncenters,
                         tree.particle_id_dtype),
                     "global_qbx_flags[i] != 0",
                     extra_args=[
@@ -857,16 +650,15 @@ class QBXFMMGeometryData(object):
         tgt_assoc = QBXTargetAssociator(self.cl_context)
 
         tgt_info = self.target_info()
-        center_info = self.center_info()
 
         from pytential.target import PointsTarget
 
         with cl.CommandQueue(self.cl_context) as queue:
             target_side_prefs = (self
-                .target_side_preferences()[center_info.ncenters:].get(queue=queue))
+                .target_side_preferences()[self.ncenters:].get(queue=queue))
 
         target_discrs_and_qbx_sides = [(
-                PointsTarget(tgt_info.targets[:, center_info.ncenters:]),
+                PointsTarget(tgt_info.targets[:, self.ncenters:]),
                 target_side_prefs.astype(np.int32))]
 
         # FIXME: try block...
@@ -878,8 +670,8 @@ class QBXFMMGeometryData(object):
 
         with cl.CommandQueue(self.cl_context) as queue:
             result = cl.array.empty(queue, tgt_info.ntargets, tree.particle_id_dtype)
-            result[:center_info.ncenters].fill(target_state.NO_QBX_NEEDED)
-            result[center_info.ncenters:] = tgt_assoc_result.target_to_center
+            result[:self.ncenters].fill(target_state.NO_QBX_NEEDED)
+            result[self.ncenters:] = tgt_assoc_result.target_to_center
 
         return result
 
@@ -891,7 +683,6 @@ class QBXFMMGeometryData(object):
         |cached|
         """
 
-        center_info = self.center_info()
         user_ttc = self.user_target_to_center()
 
         with cl.CommandQueue(self.cl_context) as queue:
@@ -917,7 +708,7 @@ class QBXFMMGeometryData(object):
             center_target_starts, targets_sorted_by_center, _ = \
                     self.code_getter.key_value_sort(queue,
                             filtered_tree_ttc, filtered_target_ids,
-                            center_info.ncenters, tree_ttc.dtype)
+                            self.ncenters, tree_ttc.dtype)
 
             logger.info("build center -> targets lookup table: done")
 
@@ -927,31 +718,6 @@ class QBXFMMGeometryData(object):
 
             return result
 
-    @memoize_method
-    def global_qbx_centers_box_target_lists(self):
-        """Build a list of targets per box consisting only of global QBX centers.
-        Returns a :class:`boxtree.tree.FilteredTargetListsInUserOrder`.
-        (I.e. no new target order is created for these targets, as we expect
-        there to be (relatively) not many of them.)
-
-        |cached|
-        """
-
-        center_info = self.center_info()
-        with cl.CommandQueue(self.cl_context) as queue:
-            logger.info("find global qbx centers box target list: start")
-
-            flags = cl.array.zeros(queue, self.tree().ntargets, np.int8)
-
-            flags[:center_info.ncenters] = self.global_qbx_flags()
-
-            from boxtree.tree import filter_target_lists_in_user_order
-            result = filter_target_lists_in_user_order(queue, self.tree(), flags)
-
-            logger.info("find global qbx centers box target list: done")
-
-            return result.with_queue(None)
-
     @memoize_method
     def non_qbx_box_target_lists(self):
         """Build a list of targets per box that don't need to bother with QBX.
@@ -973,7 +739,7 @@ class QBXFMMGeometryData(object):
 
             # 'flags' is in user order, and should be.
 
-            nqbx_centers = self.center_info().ncenters
+            nqbx_centers = self.ncenters
             flags[:nqbx_centers] = 0
 
             from boxtree.tree import filter_target_lists_in_tree_order
@@ -1012,22 +778,22 @@ class QBXFMMGeometryData(object):
 
             # {{{ draw centers and circles
 
-            center_info = self.center_info()
+            centers = self.centers()
             centers = [
-                    center_info.centers[0].get(queue),
-                    center_info.centers[1].get(queue)]
+                    centers[0].get(queue),
+                    centers[1].get(queue)]
             pt.plot(centers[0][global_flags == 0],
                     centers[1][global_flags == 0], "oc",
                     label="centers needing local qbx")
             ax = pt.gca()
-            for icenter, (cx, cy, r) in enumerate(zip(
-                    centers[0], centers[1], center_info.radii.get(queue))):
-                ax.add_artist(
-                        pt.Circle((cx, cy), r, fill=False, ls="dotted", lw=1))
-                pt.text(cx, cy,
-                        str(icenter), fontsize=8,
-                        ha="left", va="center",
-                        bbox=dict(facecolor='white', alpha=0.5, lw=0))
+            # for icenter, (cx, cy, r) in enumerate(zip(
+            #         centers[0], centers[1], center_info.radii.get(queue))):
+            #     ax.add_artist(
+            #             pt.Circle((cx, cy), r, fill=False, ls="dotted", lw=1))
+            #     pt.text(cx, cy,
+            #             str(icenter), fontsize=8,
+            #             ha="left", va="center",
+            #             bbox=dict(facecolor='white', alpha=0.5, lw=0))
 
             # }}}
 
@@ -1054,9 +820,9 @@ class QBXFMMGeometryData(object):
             tccount = 0
             checked = 0
             for tx, ty, tcenter in zip(
-                    targets[0][center_info.ncenters:],
-                    targets[1][center_info.ncenters:],
-                    ttc[center_info.ncenters:]):
+                    targets[0][self.ncenters:],
+                    targets[1][self.ncenters:],
+                    ttc[self.ncenters:]):
                 checked += 1
                 if tcenter >= 0:
                     tccount += 1
diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py
index 9ae54bb681bc96f98c31a2027b2b466868094531..8efd9a8800f614f9d54a5d30285f108c42072ba1 100644
--- a/pytential/qbx/refinement.py
+++ b/pytential/qbx/refinement.py
@@ -270,6 +270,8 @@ class RefinerWrangler(object):
         self.code_container = code_container
         self.queue = queue
 
+    # {{{ check subroutines for conditions 1-3
+
     def check_expansion_disks_undisturbed_by_sources(self,
             lpot_source, tree, peer_lists, refine_flags,
             debug, wait_for=None):
@@ -295,9 +297,9 @@ class RefinerWrangler(object):
         found_panel_to_refine.finish()
         unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
 
-        center_danger_zone_radii_by_panel = (
-                lpot_source.panel_sizes("npanels")
-                .with_queue(self.queue) / 2)
+        # FIXME: This shouldn't be by panel
+        center_danger_zone_radii_by_panel = \
+                lpot_source._expansion_radii("npanels")
 
         evt = knl(
             *unwrap_args(
@@ -354,7 +356,7 @@ class RefinerWrangler(object):
         found_panel_to_refine.finish()
 
         source_danger_zone_radii_by_panel = (
-                lpot_source.fine_panel_sizes("npanels")
+                lpot_source._fine_panel_sizes("npanels")
                 .with_queue(self.queue) / 4)
         unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
 
@@ -398,7 +400,7 @@ class RefinerWrangler(object):
             npanels_to_refine_prev = cl.array.sum(refine_flags).get()
 
         evt, out = knl(self.queue,
-                       panel_sizes=lpot_source.panel_sizes("npanels"),
+                       panel_sizes=lpot_source._panel_sizes("npanels"),
                        refine_flags=refine_flags,
                        refine_flags_updated=np.array(0),
                        kernel_length_scale=np.array(kernel_length_scale),
@@ -416,6 +418,8 @@ class RefinerWrangler(object):
 
         return (out["refine_flags_updated"].get() == 1).all()
 
+    # }}}
+
     def build_tree(self, lpot_source, use_base_fine_discr=False):
         tb = self.code_container.tree_builder()
         from pytential.qbx.utils import build_tree_with_qbx_metadata
@@ -470,6 +474,8 @@ def make_empty_refine_flags(queue, lpot_source, use_base_fine_discr=False):
     return result
 
 
+# {{{ main entry point
+
 def refine_for_global_qbx(lpot_source, code_container,
         group_factory, fine_group_factory, kernel_length_scale=None,
         # FIXME: Set debug=False once everything works.
@@ -621,6 +627,7 @@ def refine_for_global_qbx(lpot_source, code_container,
     lpot_source = lpot_source.copy(debug=debug, refined_for_global_qbx=True)
 
     if len(connections) == 0:
+        # FIXME: This is inefficient
         connection = make_same_mesh_connection(
                 lpot_source.density_discr,
                 lpot_source.density_discr)
@@ -629,5 +636,7 @@ def refine_for_global_qbx(lpot_source, code_container,
 
     return lpot_source, connection
 
+# }}}
+
 
 # vim: foldmethod=marker:filetype=pyopencl
diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py
index b5f6b153e2e340f8e4f2315be2543017a469fe2f..4255d865f753da087e1cca1e7942248d14b684b0 100644
--- a/pytential/qbx/target_assoc.py
+++ b/pytential/qbx/target_assoc.py
@@ -128,7 +128,7 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
         particle_id_t source_offset,
         particle_id_t target_offset,
         particle_id_t *sorted_target_ids,
-        coord_t *panel_sizes,
+        coord_t *tunnel_radius_by_source,
         coord_t *box_to_search_dist,
 
         /* output */
@@ -140,7 +140,7 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
             coord_t *particles_${ax},
         %endfor
     """,
-    ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""
+    ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""//CL//
         coord_vec_t tgt_coords;
         ${load_particle("INDEX_FOR_TARGET_PARTICLE(i)", "tgt_coords")}
         {
@@ -149,7 +149,7 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
             ${ball_radius} = box_to_search_dist[my_box];
         }
     """,
-    leaf_found_op=QBX_TREE_MAKO_DEFS + r"""
+    leaf_found_op=QBX_TREE_MAKO_DEFS + r"""//CL//
         for (particle_id_t source_idx = box_to_source_starts[${leaf_box_id}];
              source_idx < box_to_source_starts[${leaf_box_id} + 1];
              ++source_idx)
@@ -158,7 +158,8 @@ QBX_TARGET_MARKER = AreaQueryElementwiseTemplate(
             coord_vec_t source_coords;
             ${load_particle("INDEX_FOR_SOURCE_PARTICLE(source)", "source_coords")}
 
-            if (distance(source_coords, tgt_coords) <= panel_sizes[source] / 2)
+            if (distance(source_coords, tgt_coords)
+                    <= tunnel_radius_by_source[source])
             {
                 target_status[i] = MARKED_QBX_CENTER_PENDING;
                 *found_target_close_to_panel = 1;
@@ -177,9 +178,8 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate(
         particle_id_t center_offset,
         particle_id_t target_offset,
         particle_id_t *sorted_target_ids,
-        coord_t *panel_sizes,
+        coord_t *expansion_radii_by_center_with_stick_out,
         coord_t *box_to_search_dist,
-        coord_t stick_out_factor,
         int *target_flags,
 
         /* input/output */
@@ -222,7 +222,7 @@ QBX_CENTER_FINDER = AreaQueryElementwiseTemplate(
             coord_t my_dist_to_center = distance(tgt_coords, center_coords);
 
             if (my_dist_to_center
-                    <= (panel_sizes[center] / 2) * (1 + stick_out_factor)
+                    <= expansion_radii_by_center_with_stick_out[center]
                 && my_dist_to_center < min_dist_to_center[i])
             {
                 target_status[i] = MARKED_QBX_CENTER_FOUND;
@@ -245,7 +245,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
         particle_id_t target_offset,
         int npanels,
         particle_id_t *sorted_target_ids,
-        coord_t *panel_sizes,
+        coord_t *tunnel_radius_by_source,
         int *target_status,
         coord_t *box_to_search_dist,
 
@@ -258,7 +258,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
             coord_t *particles_${ax},
         %endfor
     """,
-    ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""
+    ball_center_and_radius_expr=QBX_TREE_C_PREAMBLE + QBX_TREE_MAKO_DEFS + r"""//CL//
         coord_vec_t target_coords;
         ${load_particle("INDEX_FOR_TARGET_PARTICLE(i)", "target_coords")}
         {
@@ -267,7 +267,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
             ${ball_radius} = box_to_search_dist[my_box];
         }
     """,
-    leaf_found_op=QBX_TREE_MAKO_DEFS + r"""
+    leaf_found_op=QBX_TREE_MAKO_DEFS + r"""//CL//
         for (particle_id_t source_idx = box_to_source_starts[${leaf_box_id}];
              source_idx < box_to_source_starts[${leaf_box_id} + 1];
              ++source_idx)
@@ -279,7 +279,7 @@ QBX_FAILED_TARGET_ASSOCIATION_REFINER = AreaQueryElementwiseTemplate(
 
             bool is_close =
                 distance(target_coords, source_coords)
-                <= panel_sizes[source] / 2;
+                <= tunnel_radius_by_source[source];
 
             if (is_close && target_status[i] == MARKED_QBX_CENTER_PENDING)
             {
@@ -398,19 +398,37 @@ class QBXTargetAssociator(object):
         # Perform a space invader query over the sources.
         source_slice = tree.user_source_ids[tree.qbx_user_source_slice]
         sources = [axis.with_queue(queue)[source_slice] for axis in tree.sources]
-        panel_sizes = lpot_source.panel_sizes("nsources").with_queue(queue)
+        tunnel_radius_by_source = \
+                lpot_source._close_target_tunnel_radius("nsources").with_queue(queue)
+
+        # Target-marking algorithm (TGTMARK):
+        #
+        # (1) Use a space invader query to tag each leaf box that intersects with the
+        # "near-source-detection tunnel" with the distance to the closest source.
+        #
+        # (2) Do an area query around all targets with the radius resulting
+        # from the space invader query, enumerate sources in that vicinity.
+        # If a source is found whose distance to the target is less than the
+        # source's tunnel radius, mark that target as pending.
+        # (or below: mark the source for refinement)
+
+        # Note that this comment is referred to below by "TGTMARK". If you
+        # remove this comment or change the algorithm here, make sure that
+        # the reference below is still accurate.
 
         box_to_search_dist, evt = self.space_invader_query(
                 queue,
                 tree,
                 sources,
-                panel_sizes / 2,
+                tunnel_radius_by_source,
                 peer_lists,
                 wait_for=wait_for)
         wait_for = [evt]
 
         logger.info("target association: marking targets close to panels")
 
+        tunnel_radius_by_source = lpot_source._close_target_tunnel_radius("nsources")
+
         evt = knl(
             *unwrap_args(
                 tree, peer_lists,
@@ -419,7 +437,7 @@ class QBXTargetAssociator(object):
                 tree.qbx_user_source_slice.start,
                 tree.qbx_user_target_slice.start,
                 tree.sorted_target_ids,
-                panel_sizes,
+                tunnel_radius_by_source,
                 box_to_search_dist,
                 target_status,
                 found_target_close_to_panel,
@@ -463,13 +481,22 @@ class QBXTargetAssociator(object):
         center_slice = \
                 tree.sorted_target_ids[tree.qbx_user_center_slice].with_queue(queue)
         centers = [axis.with_queue(queue)[center_slice] for axis in tree.sources]
-        panel_sizes = lpot_source.panel_sizes("ncenters").with_queue(queue)
+        expansion_radii_by_center = \
+                lpot_source._expansion_radii("ncenters").with_queue(queue)
+        expansion_radii_by_center_with_stick_out = \
+                expansion_radii_by_center * (1 + stick_out_factor)
+
+        # Idea:
+        #
+        # (1) Tag leaf boxes around centers with max distance to usable center.
+        # (2) Area query from targets with those radii to find closest eligible
+        # center.
 
         box_to_search_dist, evt = self.space_invader_query(
                 queue,
                 tree,
                 centers,
-                panel_sizes * ((1 + stick_out_factor) / 2),
+                expansion_radii_by_center_with_stick_out,
                 peer_lists,
                 wait_for=wait_for)
         wait_for = [evt]
@@ -490,9 +517,8 @@ class QBXTargetAssociator(object):
                 tree.qbx_user_center_slice.start,
                 tree.qbx_user_target_slice.start,
                 tree.sorted_target_ids,
-                panel_sizes,
+                expansion_radii_by_center_with_stick_out,
                 box_to_search_dist,
-                stick_out_factor,
                 target_flags,
                 target_status,
                 target_assoc.target_to_center,
@@ -535,13 +561,16 @@ class QBXTargetAssociator(object):
         # Perform a space invader query over the sources.
         source_slice = tree.user_source_ids[tree.qbx_user_source_slice]
         sources = [axis.with_queue(queue)[source_slice] for axis in tree.sources]
-        panel_sizes = lpot_source.panel_sizes("nsources").with_queue(queue)
+        tunnel_radius_by_source = \
+                lpot_source._close_target_tunnel_radius("nsources").with_queue(queue)
+
+        # See (TGTMARK) above for algorithm.
 
         box_to_search_dist, evt = self.space_invader_query(
                 queue,
                 tree,
                 sources,
-                panel_sizes / 2,
+                tunnel_radius_by_source,
                 peer_lists,
                 wait_for=wait_for)
         wait_for = [evt]
@@ -558,7 +587,7 @@ class QBXTargetAssociator(object):
                 tree.qbx_user_target_slice.start,
                 tree.nqbxpanels,
                 tree.sorted_target_ids,
-                lpot_source.panel_sizes("nsources"),
+                lpot_source._close_target_tunnel_radius("nsources"),
                 target_status,
                 box_to_search_dist,
                 refine_flags,
diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py
index 8b363250d644bfaebb65d738208ebab057c677aa..02a19e007eae86e90066224f97b0b261a02a7e6a 100644
--- a/pytential/qbx/utils.py
+++ b/pytential/qbx/utils.py
@@ -99,8 +99,8 @@ def get_interleaved_centers(queue, lpot_source):
     next to corresponding exterior centers.
     """
     knl = get_interleaver_kernel(lpot_source.density_discr.real_dtype)
-    int_centers = lpot_source.centers(-1)
-    ext_centers = lpot_source.centers(+1)
+    int_centers = get_centers_on_side(lpot_source, -1)
+    ext_centers = get_centers_on_side(lpot_source, +1)
 
     result = []
     wait_for = []
@@ -118,6 +118,160 @@ def get_interleaved_centers(queue, lpot_source):
 # }}}
 
 
+# {{{ make interleaved radii
+
+def get_interleaved_radii(queue, lpot_source):
+    """
+    Return an array of shape (dim, ncenters) in which interior centers are placed
+    next to corresponding exterior centers.
+    """
+    knl = get_interleaver_kernel(lpot_source.density_discr.real_dtype)
+    radii = lpot_source._expansion_radii("nsources")
+
+    result = cl.array.empty(queue, len(radii) * 2, radii.dtype)
+    evt, _ = knl(queue, src1=radii, src2=radii, dst=result)
+    evt.wait()
+
+    return result
+
+# }}}
+
+
+# {{{ panel sizes
+
+def panel_sizes(discr, last_dim_length):
+    if last_dim_length not in ("nsources", "ncenters", "npanels"):
+        raise ValueError(
+                "invalid value of last_dim_length: %s" % last_dim_length)
+
+    # To get the panel size this does the equivalent of (∫ 1 ds)**(1/dim).
+    # FIXME: Kernel optimizations
+
+    if last_dim_length == "nsources" or last_dim_length == "ncenters":
+        knl = lp.make_kernel(
+            "{[i,j,k]: 0<=i<nelements and 0<=j,k<nunit_nodes}",
+            "panel_sizes[i,j] = sum(k, ds[i,k])**(1/dim)",
+            name="compute_size")
+
+        def panel_size_view(discr, group_nr):
+            return discr.groups[group_nr].view
+
+    elif last_dim_length == "npanels":
+        knl = lp.make_kernel(
+            "{[i,j]: 0<=i<nelements and 0<=j<nunit_nodes}",
+            "panel_sizes[i] = sum(j, ds[i,j])**(1/dim)",
+            name="compute_size")
+        from functools import partial
+
+        def panel_size_view(discr, group_nr):
+            return partial(el_view, discr, group_nr)
+
+    else:
+        raise ValueError("unknown dim length specified")
+
+    knl = lp.fix_parameters(knl, dim=discr.dim)
+
+    with cl.CommandQueue(discr.cl_context) as queue:
+        from pytential import bind, sym
+        ds = bind(
+                discr,
+                sym.area_element(ambient_dim=discr.ambient_dim, dim=discr.dim)
+                * sym.QWeight()
+                )(queue)
+        panel_sizes = cl.array.empty(
+            queue, discr.nnodes
+            if last_dim_length in ("nsources", "ncenters")
+            else discr.mesh.nelements, discr.real_dtype)
+        for group_nr, group in enumerate(discr.groups):
+            _, (result,) = knl(queue,
+                nelements=group.nelements,
+                nunit_nodes=group.nunit_nodes,
+                ds=group.view(ds),
+                panel_sizes=panel_size_view(
+                    discr, group_nr)(panel_sizes))
+        panel_sizes.finish()
+        if last_dim_length == "ncenters":
+            from pytential.qbx.utils import get_interleaver_kernel
+            knl = get_interleaver_kernel(discr.real_dtype)
+            _, (panel_sizes,) = knl(queue, dstlen=2*discr.nnodes,
+                                    src1=panel_sizes, src2=panel_sizes)
+        return panel_sizes.with_queue(None)
+
+# }}}
+
+
+# {{{ element centers of mass
+
+def element_centers_of_mass(discr):
+    knl = lp.make_kernel(
+        """{[dim,k,i]:
+            0<=dim<ndims and
+            0<=k<nelements and
+            0<=i<nunit_nodes}""",
+        """
+            panels[dim, k] = sum(i, nodes[dim, k, i])/nunit_nodes
+            """,
+        default_offset=lp.auto, name="find_panel_centers_of_mass")
+
+    knl = lp.fix_parameters(knl, ndims=discr.ambient_dim)
+
+    knl = lp.split_iname(knl, "k", 128, inner_tag="l.0", outer_tag="g.0")
+    knl = lp.tag_inames(knl, dict(dim="ilp"))
+
+    with cl.CommandQueue(discr.cl_context) as queue:
+        mesh = discr.mesh
+        panels = cl.array.empty(queue, (mesh.ambient_dim, mesh.nelements),
+                                dtype=discr.real_dtype)
+        for group_nr, group in enumerate(discr.groups):
+            _, (result,) = knl(queue,
+                nelements=group.nelements,
+                nunit_nodes=group.nunit_nodes,
+                nodes=group.view(discr.nodes()),
+                panels=el_view(discr, group_nr, panels))
+        panels.finish()
+        panels = panels.with_queue(None)
+        return tuple(panels[d, :] for d in range(mesh.ambient_dim))
+
+# }}}
+
+
+# {{{ compute center array
+
+def get_centers_on_side(lpot_src, sign):
+    adim = lpot_src.density_discr.ambient_dim
+    dim = lpot_src.density_discr.dim
+
+    from pytential import sym, bind
+    with cl.CommandQueue(lpot_src.cl_context) as queue:
+        nodes = bind(lpot_src.density_discr, sym.nodes(adim))(queue)
+        normals = bind(lpot_src.density_discr, sym.normal(adim, dim=dim))(queue)
+        expansion_radii = lpot_src._expansion_radii("nsources").with_queue(queue)
+        return (nodes + normals * sign * expansion_radii).as_vector(np.object)
+
+# }}}
+
+
+# {{{ el_view
+
+def el_view(discr, group_nr, global_array):
+    """Return a view of *global_array* of shape
+    ``(..., discr.groups[group_nr].nelements)``
+    where *global_array* is of shape ``(..., nelements)``,
+    where *nelements* is the global (per-discretization) node count.
+    """
+
+    group = discr.groups[group_nr]
+    el_nr_base = sum(group.nelements for group in discr.groups[:group_nr])
+
+    return global_array[
+        ..., el_nr_base:el_nr_base + group.nelements] \
+        .reshape(
+            global_array.shape[:-1]
+            + (group.nelements,))
+
+# }}}
+
+
 # {{{ discr plotter
 
 def plot_discr(lpot_source, outfilename="discr.pdf"):
@@ -150,8 +304,7 @@ def plot_discr(lpot_source, outfilename="discr.pdf"):
 # }}}
 
 
-# {{{ tree creation
-
+# {{{ tree-with-metadata: data structure
 
 class TreeWithQBXMetadata(Tree):
     """A subclass of :class:`boxtree.tree.Tree`. Has all of that class's
@@ -170,41 +323,41 @@ class TreeWithQBXMetadata(Tree):
 
     .. attribute:: box_to_qbx_panel_starts
 
-    ``box_id_t [nboxes + 1]``
+        ``box_id_t [nboxes + 1]``
 
     .. attribute:: box_to_qbx_panel_lists
 
-    ``particle_id_t [*]``
+        ``particle_id_t [*]``
 
     .. rubric:: Box to QBX sources
 
     .. attribute:: box_to_qbx_source_starts
 
-    ``box_id_t [nboxes + 1]``
+        ``box_id_t [nboxes + 1]``
 
     .. attribute:: box_to_qbx_source_lists
 
-    ``particle_id_t [*]``
+        ``particle_id_t [*]``
 
     .. rubric:: Box to QBX centers
 
     .. attribute:: box_to_qbx_center_starts
 
-    ``box_id_t [nboxes + 1]``
+        ``box_id_t [nboxes + 1]``
 
     .. attribute:: box_to_qbx_center_lists
 
-    ``particle_id_t [*]``
+        ``particle_id_t [*]``
 
     .. rubric:: Box to QBX targets
 
     .. attribute:: box_to_qbx_target_starts
 
-    ``box_id_t [nboxes + 1]``
+        ``box_id_t [nboxes + 1]``
 
     .. attribute:: box_to_qbx_target_lists
 
-    ``particle_id_t [*]``
+        ``particle_id_t [*]``
 
     .. ------------------------------------------------------------------------
     .. rubric:: Panel properties
@@ -212,11 +365,11 @@ class TreeWithQBXMetadata(Tree):
 
     .. attribute:: qbx_panel_to_source_starts
 
-    ``particle_id_t [nqbxpanels + 1]``
+        ``particle_id_t [nqbxpanels + 1]``
 
     .. attribute:: qbx_panel_to_center_starts
 
-    ``particle_id_t [nqbxpanels + 1]``
+        ``particle_id_t [nqbxpanels + 1]``
 
     .. ------------------------------------------------------------------------
     .. rubric:: Particle order indices
@@ -229,6 +382,10 @@ class TreeWithQBXMetadata(Tree):
     """
     pass
 
+# }}}
+
+
+# {{{ tree-with-metadata: creation
 
 MAX_REFINE_WEIGHT = 64
 
@@ -271,9 +428,9 @@ def build_tree_with_qbx_metadata(
     centers = get_interleaved_centers(queue, lpot_source)
 
     centers_of_mass = (
-            lpot_source.panel_centers_of_mass()
+            lpot_source._panel_centers_of_mass()
             if not use_base_fine_discr
-            else lpot_source.fine_panel_centers_of_mass())
+            else lpot_source._fine_panel_centers_of_mass())
 
     targets = (tgt.nodes() for tgt in targets_list)
 
diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py
index 9b4622fde5f4bf8cc0527b05b09c38433fdf4267..2012b9928c0dad7e10947ba110c5ca65e9041faa 100644
--- a/pytential/symbolic/matrix.py
+++ b/pytential/symbolic/matrix.py
@@ -183,11 +183,13 @@ class MatrixBuilder(EvaluationMapperBase):
 
         assert target_discr is source.density_discr
 
+        from pytential.qbx.utils import get_centers_on_side
+
         assert abs(expr.qbx_forced_limit) > 0
         _, (mat,) = mat_gen(self.queue,
                 target_discr.nodes(),
                 source.fine_density_discr.nodes(),
-                source.centers(expr.qbx_forced_limit),
+                get_centers_on_side(source, expr.qbx_forced_limit),
                 **kernel_args)
 
         mat = mat.get()
diff --git a/pytential/symbolic/pde/maxwell/waveguide.py.orig b/pytential/symbolic/pde/maxwell/waveguide.py.orig
new file mode 100644
index 0000000000000000000000000000000000000000..4b5e192f698d55c36240305f3f3d7bc810669417
--- /dev/null
+++ b/pytential/symbolic/pde/maxwell/waveguide.py.orig
@@ -0,0 +1,989 @@
+from __future__ import division, absolute_import, print_function
+
+__copyright__ = "Copyright (C) 2013-2016 Andreas Kloeckner"
+
+__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.
+"""
+
+
+__doc__ = """
+Second-Kind Waveguide
+^^^^^^^^^^^^^^^^^^^^^
+
+.. autoclass:: SecondKindInfZMuellerOperator
+
+2D Dielectric (old-style)
+^^^^^^^^^^^^^^^^^^^^^^^^^
+
+.. autoclass:: DielectricSRep2DBoundaryOperator
+.. autoclass:: DielectricSDRep2DBoundaryOperator
+"""
+
+import numpy as np
+from collections import namedtuple
+from six.moves import range
+
+from pytential import sym
+from pytools import memoize_method
+from pytential.symbolic.pde.scalar import L2WeightedPDEOperator
+
+
+# {{{ second-kind infinite-z Mueller operator
+
+class SecondKindInfZMuellerOperator(L2WeightedPDEOperator):
+    """
+    Second-kind IE representation by
+    `Lai and Jiang <https://arxiv.org/abs/1512.01117>`_.
+    """
+
+    def __init__(self, domain_n_exprs, ne,
+            interfaces, use_l2_weighting=None):
+        """
+        :attr interfaces: a tuple of tuples
+            ``(outer_domain, inner_domain, interface_id)``,
+            where *outer_domain* and *inner_domain* are indices into
+            *domain_k_names*,
+            and *interface_id* is a symbolic name for the discretization of the
+            interface. 'outer' designates the side of the interface to which
+            the normal points.
+        :attr domain_n_exprs: a tuple of variable names of the Helmholtz
+            parameter *k*, to be used inside each part of the source geometry.
+            May also be a tuple of strings, which will be transformed into
+            variable references of the corresponding names.
+        :attr beta: A symbolic expression for the wave number in the :math:`z`
+            direction. May be a string, which will be interpreted as a variable
+            name.
+        """
+
+        self.interfaces = interfaces
+
+        ne = sym.var(ne)
+        self.ne = sym.cse(ne, "ne")
+
+        self.domain_n_exprs = [
+                sym.var(n_expr)
+                for idom, n_expr in enumerate(domain_n_exprs)]
+        del domain_n_exprs
+
+        import pymbolic.primitives as p
+
+        def upper_half_square_root(x):
+            return p.If(
+                    p.Comparison(
+                        (x**0.5).a.imag,
+                        "<", 0),
+                    1j*(-x)**0.5,
+                    x**0.5)
+
+        self.domain_K_exprs = [
+                sym.cse(
+                    upper_half_square_root(n_expr**2-ne**2),
+                    "K%d" % i)
+                for i, n_expr in enumerate(self.domain_n_exprs)]
+
+        from sumpy.kernel import HelmholtzKernel
+        self.kernel = HelmholtzKernel(2, allow_evanescent=True)
+
+    def make_unknown(self, name):
+        return sym.make_sym_vector(name, len(self.unknown_index))
+
+    @property
+    @memoize_method
+    def unknown_index(self):
+        result = {}
+        for i in range(len(self.interfaces)):
+            for current in ["j", "m"]:
+                for cur_dir in ["z", "t"]:
+                    result[(current + cur_dir, i)] = len(result)
+
+        return result
+
+    def tangent(self, where):
+        return sym.cse(
+                (sym.pseudoscalar(2, 1, where)
+                / sym.area_element(2, 1, where)),
+                "tangent")
+
+    def S(self, dom_idx, density, qbx_forced_limit=+1):  # noqa
+        return sym.S(
+                self.kernel, density,
+                k=self.domain_K_exprs[dom_idx],
+                qbx_forced_limit=qbx_forced_limit)
+
+    def D(self, dom_idx, density, qbx_forced_limit="avg"):  # noqa
+        return sym.D(
+                self.kernel, density,
+                k=self.domain_K_exprs[dom_idx],
+                qbx_forced_limit=qbx_forced_limit)
+
+    def T(self, where, dom_idx, density):  # noqa
+        return sym.int_g_dsource(
+                2,
+                self.tangent(where),
+                self.kernel,
+                density,
+                k=self.domain_K_exprs[dom_idx],
+                # ????
+                qbx_forced_limit=1).xproject(0)
+
+    def representation(self, domain_idx, unknown):
+        """"Return a tuple of vectors [Ex, Ey, Ez] and [Hx, Hy, Hz]
+        representing the solution to Maxwell's equation on domain
+        *domain_idx*.
+        """
+        result = np.zeros(4, dtype=object)
+
+        unk_idx = self.unknown_index
+        k_v = self.k_vacuum
+
+        for i in range(len(self.interfaces)):
+            dom0i, dom1i, where = self.interfaces[i]
+
+            if domain_idx == dom0i:
+                domi = dom0i
+            elif domain_idx == dom1i:
+                domi = dom1i
+            else:
+                # Interface does not border the requested domain
+                continue
+
+            beta = self.ne*k_v
+            k = sym.cse(k_v * self.domain_n_exprs[domi],
+                    "k%d" % domi)
+
+            jt = unknown[unk_idx["jz", i]]
+            jz = unknown[unk_idx["jt", i]]
+            mt = unknown[unk_idx["mz", i]]
+            mz = unknown[unk_idx["mt", i]]
+
+            def diff_axis(iaxis, operand):
+                v = np.array(2, dtype=np.float64)
+                v[iaxis] = 1
+                d = sym.Derivative()
+                return d.resolve(
+                        (sym.MultiVector(v).scalar_product(d.dnabla(2)))
+                        * d(operand))
+
+            from functools import partial
+            dx = partial(diff_axis, 1)
+            dy = partial(diff_axis, 1)
+
+            tangent = self.tangent(where)
+            tau1 = tangent[0]
+            tau2 = tangent[1]
+
+            S = self.S  # noqa
+            D = self.D  # noqa
+            T = self.T  # noqa
+
+            # Ex
+            result[0] += (
+                    -1/(1j*k_v) * dx(T(where, domi, jt))
+                    + beta/k_v * dx(S(domi, jz))
+                    + k**2/(1j*k_v) * S(domi, jt * tau1)
+                    - dy(S(domi, mz))
+                    + 1j*beta * S(domi, mt*tau2)
+                    )
+
+            # Ey
+            result[1] += (
+                    - 1/(1j*k_v) * dy(T(where, domi, jt))
+                    + beta/k_v * dy(S(domi, jz))
+                    + k**2/(1j*k_v) * S(domi, jt * tau2)
+                    + dx(S(domi, mz))
+                    - 1j*beta*S(domi, mt * tau1)
+                    )
+
+            # Ez
+            result[2] += (
+                    - beta/k_v * T(where, domi, jt)
+                    + (k**2 - beta**2)/(1j*k_v) * S(domi, jz)
+                    + D(domi, mt)
+                    )
+
+            # Hx
+            result[3] += (
+                    1/(1j*k_v) * dx(T(where, domi, mt))
+                    - beta/k_v * dx(S(domi, mz))
+                    - k**2/(1j*k_v) * S(domi, mt*tau1)
+                    - k**2/k_v**2 * dy(S(domi, jz))
+                    + 1j * beta * k**2/k_v**2 * S(domi, jt*tau2)
+                    )
+
+            # Hy
+            result[4] += (
+                    1/(1j*k_v) * dy(T(where, domi, mt))
+                    - beta/k_v * dy(S(domi, mz))
+                    - k**2/(1j*k_v) * S(domi, mt * tau2)
+                    + k**2/k_v**2 * dx(S(domi, jz))
+                    - 1j*beta * k**2/k_v**2 * S(domi, jt*tau1)
+                    )
+
+            # Hz
+            result[5] += (
+                    beta/k_v * T(where, domi, mt)
+                    - (k**2 - beta**2)/(1j*k_v) * S(domi, mz)
+                    + k**2/k_v**2 * D(domi, jt)
+                    )
+
+        return result
+
+    def operator(self, unknown):
+        result = np.zeros(4*len(self.interfaces), dtype=object)
+
+        unk_idx = self.unknown_index
+
+        for i in range(len(self.interfaces)):
+            idx_jt = unk_idx["jz", i]
+            idx_jz = unk_idx["jt", i]
+            idx_mt = unk_idx["mz", i]
+            idx_mz = unk_idx["mt", i]
+
+            phi1 = unknown[idx_jt]
+            phi2 = unknown[idx_jz]
+            phi3 = unknown[idx_mt]
+            phi4 = unknown[idx_mz]
+
+            ne = self.ne
+
+            dom0i, dom1i, where = self.interfaces[i]
+
+            tangent = self.tangent(where)
+            normal = sym.cse(
+                    sym.normal(2, 1, where),
+                    "normal")
+
+            S = self.S  # noqa
+            D = self.D  # noqa
+            T = self.T  # noqa
+
+            def Tt(where, dom, density):  # noqa
+                return sym.tangential_derivative(
+                        2, T(where, dom, density)).xproject(0)
+
+            def Sn(dom, density):  # noqa
+                return sym.normal_derivative(
+                        2,
+                        S(dom, density,
+                            qbx_forced_limit="avg"))
+
+            def St(dom, density):  # noqa
+                return sym.tangential_derivative(2, S(dom, density)).xproject(0)
+
+            n0 = self.domain_n_exprs[dom0i]
+            n1 = self.domain_n_exprs[dom1i]
+
+            a11 = sym.cse(n0**2 * D(dom0i, phi1) - n1**2 * D(dom1i, phi1), "a11")
+            a22 = sym.cse(-n0**2 * Sn(dom0i, phi2) + n1**2 * Sn(dom1i, phi2), "a22")
+            a33 = sym.cse(D(dom0i, phi3)-D(dom1i, phi3), "a33")
+            a44 = sym.cse(-Sn(dom0i, phi4) + Sn(dom1i, phi4), "a44")
+
+            a21 = sym.cse(-1j * ne * (
+                    n0**2 * tangent.scalar_product(
+                        S(dom0i, normal * phi1))
+                    - n1**2 * tangent.scalar_product(
+                        S(dom1i, normal * phi1))), "a21")
+
+            a43 = sym.cse(-1j * ne * (
+                    tangent.scalar_product(
+                        S(dom0i, normal * phi3))
+                    - tangent.scalar_product(
+                        S(dom1i, normal * phi3))), "a43")
+
+            a13 = +1*sym.cse(
+                    ne*(T(where, dom0i, phi3) - T(where, dom1i, phi3)), "a13")
+            a31 = -1*sym.cse(
+                    ne*(T(where, dom0i, phi1) - T(where, dom1i, phi1)), "a31")
+
+            a24 = +1*sym.cse(ne*(St(dom0i, phi4) - St(dom1i, phi4)), "a24")
+            a42 = -1*sym.cse(ne*(St(dom0i, phi2) - St(dom1i, phi2)), "a42")
+
+            a14 = sym.cse(1j*(
+                    (n0**2 - ne**2) * S(dom0i, phi4)
+                    - (n1**2 - ne**2) * S(dom1i, phi4)
+                    ), "a14")
+            a32 = -sym.cse(1j*(
+                    (n0**2 - ne**2) * S(dom0i, phi2)
+                    - (n1**2 - ne**2) * S(dom1i, phi2)
+                    ), "a32")
+
+            def a23_expr(phi):
+                return (
+                        1j * (Tt(where, dom0i, phi) - Tt(where, dom1i, phi))
+                        - 1j * (
+                            n0**2 * tangent.scalar_product(
+<<<<<<< HEAD
+                                S(dom0i, normal * phi))
+                            - n1**2 * tangent.scalar_product(
+                                S(dom1i, normal * phi))))
+=======
+                                S(0, tangent * phi))
+                            - n1**2 * tangent.scalar_product(
+                                S(1, tangent * phi))))
+>>>>>>> 9af8f2bde0070ee8488c6bc58c8ece5114004a71
+
+            a23 = +1*sym.cse(a23_expr(phi3), "a23")
+            a41 = -1*sym.cse(a23_expr(phi1), "a41")
+
+            d1 = (n0**2 + n1**2)/2 * phi1
+            d2 = (n0**2 + n1**2)/2 * phi2
+            d3 = phi3
+            d4 = phi4
+
+            result[idx_jt] += d1 + a11 + 000 + a13 + a14
+            result[idx_jz] += d2 + a21 + a22 + a23 + a24
+            result[idx_mt] += d3 + a31 + a32 + a33 + 0
+            result[idx_mz] += d4 + a41 + a42 + a43 + a44
+
+            # TODO: L2 weighting
+            # TODO: Add representation contributions to other boundaries
+            # abutting the domain
+            return result
+
+# }}}
+
+
+# {{{ old-style waveguide
+
+class Dielectric2DBoundaryOperatorBase(L2WeightedPDEOperator):
+    r"""
+    Solves the following system of BVPs on :math:`\mathbb{R}^2`, in which
+    a disjoint family of domains :math:`\Omega_i` is embedded:
+
+    .. math::
+
+        \triangle E + (k_i^2-\beta^2) E = 0\quad \text{on $\Omega_i$}\\
+        \triangle H + (k_i^2-\beta^2) H = 0\quad \text{on $\Omega_i$}\\
+        [H] = 0 \text{ on $\partial \Omega_i$},\quad
+        [E] = 0 \text{ on $\partial \Omega_i$}\\
+        \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}H\right] = 0
+        \quad\text{on $\partial \Omega_i$},\quad\\
+        \left[ \frac{k_0}{k^2-\beta^2} \partial_{\hat n}E\right] = 0
+        \quad\text{on $\partial \Omega_i$}
+
+    :math:`E` and :math:`H` are assumed to be of the form
+
+    .. math::
+
+        E(x,y,z,t)=E(x,y)e^{i(\beta z-\omega t)
+        H(x,y,z,t)=H(x,y)e^{i(\beta z-\omega t)
+
+    where :math:`[\cdot]` denotes the jump across an interface, and :math:`k`
+    (without an index) denotes the value of :math:`k` on either side of the
+    interface, for the purpose of computing the jump. :math:`\hat n` denotes
+    the unit normal of the interface.
+
+    .. automethod:: make_unknown
+    .. automethod:: representation_outer
+    .. automethod:: representation_inner
+    .. automethod:: operator
+    """
+
+    field_kind_e = 0
+    field_kind_h = 1
+    field_kinds = [field_kind_e, field_kind_h]
+
+    side_in = 0
+    side_out = 1
+    sides = [side_in, side_out]
+    side_to_sign = {
+            side_in: -1,
+            side_out: 1,
+            }
+
+    dir_none = 0
+    dir_normal = 1
+    dir_tangential = 2
+
+    BCTermDescriptor = namedtuple("BCDescriptor",
+            "i_interface direction field_kind coeff_inner coeff_outer".split())
+
+    # {{{ constructor
+
+    def __init__(self, mode, k_vacuum, domain_k_exprs, beta,
+            interfaces, use_l2_weighting=None):
+        """
+        :attr mode: one of 'te', 'tm', 'tem'
+        :attr k_vacuum: A symbolic expression for the wave number in vacuum.
+            May be a string, which will be interpreted as a variable name.
+        :attr interfaces: a tuple of tuples
+            ``(outer_domain, inner_domain, interface_id)``,
+            where *outer_domain* and *inner_domain* are indices into
+            *domain_k_names*,
+            and *interface_id* is a symbolic name for the discretization of the
+            interface. 'outer' designates the side of the interface to which
+            the normal points.
+        :attr domain_k_exprs: a tuple of variable names of the Helmholtz
+            parameter *k*, to be used inside each part of the source geometry.
+            May also be a tuple of strings, which will be transformed into
+            variable references of the corresponding names.
+        :attr beta: A symbolic expression for the wave number in the :math:`z`
+            direction. May be a string, which will be interpreted as a variable
+            name.
+        """
+
+        if use_l2_weighting is None:
+            use_l2_weighting = False
+
+        super(Dielectric2DBoundaryOperatorBase, self).__init__(
+                use_l2_weighting=use_l2_weighting)
+
+        if mode == "te":
+            self.ez_enabled = False
+            self.hz_enabled = True
+        elif mode == "tm":
+            self.ez_enabled = True
+            self.hz_enabled = False
+        elif mode == "tem":
+            self.ez_enabled = True
+            self.hz_enabled = True
+        else:
+            raise ValueError("invalid mode '%s'" % mode)
+
+        self.interfaces = interfaces
+
+        fk_e = self.field_kind_e
+        fk_h = self.field_kind_h
+
+        dir_none = self.dir_none
+        dir_normal = self.dir_normal
+        dir_tangential = self.dir_tangential
+
+        if isinstance(beta, str):
+            beta = sym.var(beta)
+        beta = sym.cse(beta, "beta")
+
+        if isinstance(k_vacuum, str):
+            k_vacuum = sym.var(k_vacuum)
+        k_vacuum = sym.cse(k_vacuum, "k_vac")
+
+        self.domain_k_exprs = [
+                sym.var(k_expr)
+                if isinstance(k_expr, str)
+                else sym.cse(k_expr, "k%d" % idom)
+                for idom, k_expr in enumerate(domain_k_exprs)]
+        del domain_k_exprs
+
+        # Note the case of k/K!
+        # "K" is the 2D Helmholtz parameter.
+        # "k" is the 3D Helmholtz parameter.
+
+        self.domain_K_exprs = [
+                sym.cse((k_expr**2-beta**2)**0.5, "K%d" % i)
+                for i, k_expr in enumerate(self.domain_k_exprs)]
+
+        from sumpy.kernel import HelmholtzKernel
+        self.kernel = HelmholtzKernel(2, allow_evanescent=True)
+
+        # {{{ build bc list
+
+        # list of tuples, where each tuple consists of BCTermDescriptor instances
+
+        all_bcs = []
+        for i_interface, (outer_domain, inner_domain, _) in (
+                enumerate(self.interfaces)):
+            k_outer = self.domain_k_exprs[outer_domain]
+            k_inner = self.domain_k_exprs[inner_domain]
+
+            all_bcs += [
+                    (  # [E] = 0
+                        self.BCTermDescriptor(
+                            i_interface=i_interface,
+                            direction=dir_none,
+                            field_kind=fk_e,
+                            coeff_outer=1,
+                            coeff_inner=-1),
+                        ),
+                    (  # [H] = 0
+                        self.BCTermDescriptor(
+                            i_interface=i_interface,
+                            direction=dir_none,
+                            field_kind=fk_h,
+                            coeff_outer=1,
+                            coeff_inner=-1),
+                        ),
+                    (
+                        self.BCTermDescriptor(
+                            i_interface=i_interface,
+                            direction=dir_tangential,
+                            field_kind=fk_e,
+                            coeff_outer=beta/(k_outer**2-beta**2),
+                            coeff_inner=-beta/(k_inner**2-beta**2)),
+                        self.BCTermDescriptor(
+                            i_interface=i_interface,
+                            direction=dir_normal,
+                            field_kind=fk_h,
+                            coeff_outer=sym.cse(-k_vacuum/(k_outer**2-beta**2)),
+                            coeff_inner=sym.cse(k_vacuum/(k_inner**2-beta**2))),
+                        ),
+                    (
+                        self.BCTermDescriptor(
+                            i_interface=i_interface,
+                            direction=dir_tangential,
+                            field_kind=fk_h,
+                            coeff_outer=beta/(k_outer**2-beta**2),
+                            coeff_inner=-beta/(k_inner**2-beta**2)),
+                        self.BCTermDescriptor(
+                            i_interface=i_interface,
+                            direction=dir_normal,
+                            field_kind=fk_e,
+                            coeff_outer=sym.cse(
+                                (k_outer**2/k_vacuum)/(k_outer**2-beta**2)),
+                            coeff_inner=sym.cse(
+                                -(k_inner**2/k_vacuum)
+                                / (k_inner**2-beta**2)))
+                        ),
+                    ]
+
+            del k_outer
+            del k_inner
+
+        self.bcs = []
+        for bc in all_bcs:
+            any_significant_e = any(
+                    term.field_kind == fk_e
+                    and term.direction in [dir_normal, dir_none]
+                    for term in bc)
+            any_significant_h = any(
+                    term.field_kind == fk_h
+                    and term.direction in [dir_normal, dir_none]
+                    for term in bc)
+            is_necessary = (
+                    (self.ez_enabled and any_significant_e)
+                    or
+                    (self.hz_enabled and any_significant_h))
+
+            # Only keep tangential modes for TEM. Otherwise,
+            # no jump in H already implies jump condition on
+            # tangential derivative.
+            is_tem = self.ez_enabled and self.hz_enabled
+            terms = tuple(
+                    term
+                    for term in bc
+                    if term.direction != dir_tangential
+                    or is_tem)
+
+            if is_necessary:
+                self.bcs.append(terms)
+
+        assert (len(all_bcs)
+                * (int(self.ez_enabled) + int(self.hz_enabled)) // 2
+                == len(self.bcs))
+
+        # }}}
+
+    # }}}
+
+    def is_field_present(self, field_kind):
+        return (
+                (field_kind == self.field_kind_e and self.ez_enabled)
+                or
+                (field_kind == self.field_kind_h and self.hz_enabled))
+
+    def make_unknown(self, name):
+        num_densities = (
+                2
+                * (int(self.ez_enabled) + int(self.hz_enabled))
+                * len(self.interfaces))
+
+        assert num_densities == len(self.bcs)
+
+        return sym.make_sym_vector(name, num_densities)
+
+    def bc_term_to_operator_contrib(self, term, side, raw_potential_op,
+            density, discrete):
+        potential_op = raw_potential_op
+
+        side_sign = self.side_to_sign[side]
+
+        domain_outer, domain_inner, interface_id = \
+                self.interfaces[term.i_interface]
+        if side == self.side_in:
+            K_expr = self.domain_K_exprs[domain_inner]  # noqa
+            bc_coeff = term.coeff_inner
+        elif side == self.side_out:
+            K_expr = self.domain_K_exprs[domain_outer]  # noqa
+            bc_coeff = term.coeff_outer
+        else:
+            raise ValueError("invalid value of 'side'")
+
+        potential_op = potential_op(
+                self.kernel, density, source=interface_id,
+                k=K_expr)
+
+        if term.direction == self.dir_none:
+            if raw_potential_op is sym.S:
+                jump_term = 0
+            elif raw_potential_op is sym.D:
+                jump_term = (side_sign*0.5) * discrete
+            else:
+                assert False, raw_potential_op
+        elif term.direction == self.dir_normal:
+            potential_op = sym.normal_derivative(
+                    potential_op, interface_id)
+
+            if raw_potential_op is sym.S:
+                # S'
+                jump_term = (-side_sign*0.5) * discrete
+            elif raw_potential_op is sym.D:
+                jump_term = 0
+            else:
+                assert False, raw_potential_op
+
+        elif term.direction == self.dir_tangential:
+            potential_op = sym.tangential_derivative(
+                    raw_potential_op(
+                        self.kernel, density, source=interface_id,
+                        k=K_expr, qbx_forced_limit=side_sign),
+                    interface_id).a.as_scalar()
+
+            # Some of these may have jumps, but QBX does the dirty
+            # work here by directly computing the limit.
+            jump_term = 0
+
+        else:
+            raise ValueError("invalid direction")
+
+        potential_op = (
+                jump_term
+                + self.get_sqrt_weight(interface_id)*potential_op)
+
+        del jump_term
+
+        contrib = bc_coeff * potential_op
+
+        if (raw_potential_op is sym.D
+                and term.direction == self.dir_normal):
+            # FIXME The hypersingular part should perhaps be
+            # treated specially to avoid cancellation.
+            pass
+
+        return contrib
+
+
+# {{{ single-layer representation
+
+class DielectricSRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase):
+    def _structured_unknown(self, unknown, with_l2_weights):
+        """
+        :arg with_l2_weights: If True, return the 'bare' unknowns
+            that do not have the :math:`L^2` weights divided out.
+            Note: Those unknowns should *not* be interpreted as
+            point values of a density.
+        :returns: an array of unknowns, with the following index axes:
+            ``[side, field_kind, i_interface]``, where
+            ``side`` is o for the outside part and i for the interior part,
+            ``field_kind`` is 0 for the E-field and 1 for the H-field part,
+            ``i_interface`` is the number of the enclosed domain, starting from 0.
+        """
+        result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object)
+
+        i_unknown = 0
+        for side in self.sides:
+            for field_kind in self.field_kinds:
+                for i_interface in range(len(self.interfaces)):
+
+                    if self.is_field_present(field_kind):
+                        dens = unknown[i_unknown]
+                        i_unknown += 1
+                    else:
+                        dens = 0
+
+                    _, _, interface_id = self.interfaces[i_interface]
+
+                    if not with_l2_weights:
+                        dens = sym.cse(
+                                dens/self.get_sqrt_weight(interface_id),
+                                "dens_{side}_{field}_{dom}".format(
+                                    side={
+                                        self.side_out: "o",
+                                        self.side_in: "i"}
+                                    [side],
+                                    field={
+                                        self.field_kind_e: "E",
+                                        self.field_kind_h: "H"
+                                        }
+                                    [field_kind],
+                                    dom=i_interface))
+
+                    result[side, field_kind, i_interface] = dens
+
+        assert i_unknown == len(unknown)
+        return result
+
+    def representation(self, unknown, i_domain):
+        """
+        :return: a symbolic expression for the representation of the PDE solution
+            in domain number *i_domain*.
+        """
+        unk = self._structured_unknown(unknown, with_l2_weights=False)
+
+        result = []
+
+        for field_kind in self.field_kinds:
+            if not self.is_field_present(field_kind):
+                continue
+
+            field_result = 0
+            for i_interface, (i_domain_outer, i_domain_inner, interface_id) in (
+                    enumerate(self.interfaces)):
+                if i_domain_outer == i_domain:
+                    side = self.side_out
+                elif i_domain_inner == i_domain:
+                    side = self.side_in
+                else:
+                    continue
+
+                my_unk = unk[side, field_kind, i_interface]
+                if my_unk:
+                    field_result += sym.S(
+                            self.kernel,
+                            my_unk,
+                            source=interface_id,
+                            k=self.domain_K_exprs[i_domain])
+
+            result.append(field_result)
+
+        from pytools.obj_array import make_obj_array
+        return make_obj_array(result)
+
+    def operator(self, unknown):
+        density_unk = self._structured_unknown(unknown, with_l2_weights=False)
+        discrete_unk = self._structured_unknown(unknown, with_l2_weights=True)
+
+        result = []
+        for bc in self.bcs:
+            op = 0
+
+            for side in self.sides:
+                for term in bc:
+                    unk_index = (side, term.field_kind, term.i_interface)
+                    density = density_unk[unk_index]
+                    discrete = discrete_unk[unk_index]
+
+                    op += self.bc_term_to_operator_contrib(
+                            term, side, sym.S, density, discrete)
+
+            result.append(op)
+
+        return np.array(result, dtype=np.object)
+
+# }}}
+
+
+# {{{ single + double layer representation
+
+class DielectricSDRep2DBoundaryOperator(Dielectric2DBoundaryOperatorBase):
+    pot_kind_S = 0
+    pot_kind_D = 1
+    pot_kinds = [pot_kind_S, pot_kind_D]
+    potential_ops = {
+            pot_kind_S: sym.S,
+            pot_kind_D: sym.D,
+            }
+
+    def __init__(self, mode, k_vacuum, domain_k_exprs, beta,
+            interfaces, use_l2_weighting=None):
+
+        super(DielectricSDRep2DBoundaryOperator, self).__init__(
+                mode, k_vacuum, domain_k_exprs, beta,
+                interfaces, use_l2_weighting=use_l2_weighting)
+
+        side_in = self.side_in
+        side_out = self.side_out
+
+        def find_normal_derivative_bc_coeff(field_kind, i_interface, side):
+            result = 0
+            for bc in self.bcs:
+                for term in bc:
+                    if term.field_kind != field_kind:
+                        continue
+                    if term.i_interface != i_interface:
+                        continue
+                    if term.direction != self.dir_normal:
+                        continue
+
+                    if side == side_in:
+                        result += term.coeff_inner
+                    elif side == side_out:
+                        result += term.coeff_outer
+                    else:
+                        raise ValueError("invalid side")
+
+            return result
+
+        self.density_coeffs = np.zeros(
+                (len(self.pot_kinds), len(self.field_kinds),
+                    len(self.interfaces), len(self.sides)),
+                dtype=np.object)
+        for field_kind in self.field_kinds:
+            for i_interface in range(len(self.interfaces)):
+                self.density_coeffs[
+                        self.pot_kind_S, field_kind, i_interface, side_in] = 1
+                self.density_coeffs[
+                        self.pot_kind_S, field_kind, i_interface, side_out] = 1
+
+                # These need to satisfy
+                #
+                # [dens_coeff_D * bc_coeff * dn D]
+                # = dens_coeff_D_out * bc_coeff_out * (dn D)
+                #   + dens_coeff_D_in * bc_coeff_in * dn D
+                # = 0
+                #
+                # (because dn D is hypersingular, which we'd like to cancel out)
+                #
+                # NB: bc_coeff_{in,out} already contain the signs to realize
+                # the subtraction for the jump. (So the "+" above is as it
+                # should be.)
+
+                dens_coeff_D_in = find_normal_derivative_bc_coeff(  # noqa
+                        field_kind, i_interface, side_out)
+                dens_coeff_D_out = - find_normal_derivative_bc_coeff(  # noqa
+                        field_kind, i_interface, side_in)
+
+                self.density_coeffs[
+                        self.pot_kind_D, field_kind, i_interface, side_in] \
+                                = dens_coeff_D_in
+                self.density_coeffs[
+                        self.pot_kind_D, field_kind, i_interface, side_out] \
+                                = dens_coeff_D_out
+
+    def _structured_unknown(self, unknown, with_l2_weights):
+        """
+        :arg with_l2_weights: If True, return the 'bare' unknowns
+            that do not have the :math:`L^2` weights divided out.
+            Note: Those unknowns should *not* be interpreted as
+            point values of a density.
+        :returns: an array of unknowns, with the following index axes:
+            ``[pot_kind, field_kind, i_interface]``, where
+            ``pot_kind`` is 0 for the single-layer part and 1 for the double-layer
+            part,
+            ``field_kind`` is 0 for the E-field and 1 for the H-field part,
+            ``i_interface`` is the number of the enclosed domain, starting from 0.
+        """
+        result = np.zeros((2, 2, len(self.interfaces)), dtype=np.object)
+
+        i_unknown = 0
+        for pot_kind in self.pot_kinds:
+            for field_kind in self.field_kinds:
+                for i_interface in range(len(self.interfaces)):
+
+                    if self.is_field_present(field_kind):
+                        dens = unknown[i_unknown]
+                        i_unknown += 1
+                    else:
+                        dens = 0
+
+                    _, _, interface_id = self.interfaces[i_interface]
+
+                    if not with_l2_weights:
+                        dens = sym.cse(
+                                dens/self.get_sqrt_weight(interface_id),
+                                "dens_{pot}_{field}_{intf}".format(
+                                    pot={0: "S", 1: "D"}[pot_kind],
+                                    field={
+                                        self.field_kind_e: "E",
+                                        self.field_kind_h: "H"
+                                        }
+                                    [field_kind],
+                                    intf=i_interface))
+
+                    result[pot_kind, field_kind, i_interface] = dens
+
+        assert i_unknown == len(unknown)
+        return result
+
+    def representation(self, unknown, i_domain):
+        """
+        :return: a symbolic expression for the representation of the PDE solution
+            in domain number *i_domain*.
+        """
+        unk = self._structured_unknown(unknown, with_l2_weights=False)
+
+        result = []
+
+        for field_kind in self.field_kinds:
+            if not self.is_field_present(field_kind):
+                continue
+
+            field_result = 0
+            for pot_kind in self.pot_kinds:
+                for i_interface, (i_domain_outer, i_domain_inner, interface_id) in (
+                        enumerate(self.interfaces)):
+                    if i_domain_outer == i_domain:
+                        side = self.side_out
+                    elif i_domain_inner == i_domain:
+                        side = self.side_in
+                    else:
+                        continue
+
+                    my_unk = unk[pot_kind, field_kind, i_interface]
+                    if my_unk:
+                        field_result += (
+                                self.density_coeffs[
+                                    pot_kind, field_kind, i_interface, side]
+                                * self.potential_ops[pot_kind](
+                                    self.kernel,
+                                    my_unk,
+                                    source=interface_id,
+                                    k=self.domain_K_exprs[i_domain]
+                                    ))
+
+            result.append(field_result)
+
+        from pytools.obj_array import make_obj_array
+        return make_obj_array(result)
+
+    def operator(self, unknown):
+        density_unk = self._structured_unknown(unknown, with_l2_weights=False)
+        discrete_unk = self._structured_unknown(unknown, with_l2_weights=True)
+
+        result = []
+        for bc in self.bcs:
+            op = 0
+
+            for pot_kind in self.pot_kinds:
+                for term in bc:
+
+                    for side in self.sides:
+                        raw_potential_op = \
+                                self.potential_ops[pot_kind]
+
+                        unk_index = (pot_kind, term.field_kind, term.i_interface)
+                        density = density_unk[unk_index]
+                        discrete = discrete_unk[unk_index]
+
+                        op += (
+                                self.density_coeffs[
+                                    pot_kind, term.field_kind, term.i_interface,
+                                    side]
+                                * self.bc_term_to_operator_contrib(
+                                    term, side, raw_potential_op, density, discrete)
+                                )
+
+            result.append(op)
+
+        return np.array(result, dtype=np.object)
+
+# }}}
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py
index 4013b8f31081d6e0d4deda9a846fd4bf0e0af552..7326efc13d212be30c83adecedc7a3e579722605 100644
--- a/test/test_global_qbx.py
+++ b/test/test_global_qbx.py
@@ -107,20 +107,21 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
             lpot_source, RefinerCodeContainer(cl_ctx),
             factory, fine_factory, **refiner_extra_kwargs)
 
+    from pytential.qbx.utils import get_centers_on_side
+
     discr_nodes = lpot_source.density_discr.nodes().get(queue)
     fine_discr_nodes = lpot_source.fine_density_discr.nodes().get(queue)
-    int_centers = lpot_source.centers(-1)
+    int_centers = get_centers_on_side(lpot_source, -1)
     int_centers = np.array([axis.get(queue) for axis in int_centers])
-    ext_centers = lpot_source.centers(+1)
+    ext_centers = get_centers_on_side(lpot_source, +1)
     ext_centers = np.array([axis.get(queue) for axis in ext_centers])
-    panel_sizes = lpot_source.panel_sizes("npanels").get(queue)
-    fine_panel_sizes = lpot_source.fine_panel_sizes("npanels").get(queue)
+    expansion_radii = lpot_source._expansion_radii("npanels").get(queue)
+    panel_sizes = lpot_source._panel_sizes("npanels").get(queue)
+    fine_panel_sizes = lpot_source._fine_panel_sizes("npanels").get(queue)
 
     # {{{ check if satisfying criteria
 
     def check_disk_undisturbed_by_sources(centers_panel, sources_panel):
-        h = panel_sizes[centers_panel.element_nr]
-
         if centers_panel.element_nr == sources_panel.element_nr:
             # Same panel
             return
@@ -143,8 +144,9 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
         # A center cannot be closer to another panel than to its originating
         # panel.
 
-        assert dist >= h / 2, \
-                (dist, h, centers_panel.element_nr, sources_panel.element_nr)
+        rad = expansion_radii[centers_panel.element_nr]
+        assert dist >= rad, \
+                (dist, rad, centers_panel.element_nr, sources_panel.element_nr)
 
     def check_sufficient_quadrature_resolution(centers_panel, sources_panel):
         h = fine_panel_sizes[sources_panel.element_nr]
@@ -232,8 +234,10 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
     lpot_source, conn = QBXLayerPotentialSource(discr, order).with_refinement()
     del discr
 
-    int_centers = lpot_source.centers(-1)
-    ext_centers = lpot_source.centers(+1)
+    from pytential.qbx.utils import get_centers_on_side
+
+    int_centers = get_centers_on_side(lpot_source, -1)
+    ext_centers = get_centers_on_side(lpot_source, +1)
 
     # }}}
 
@@ -243,7 +247,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
     rng = PhiloxGenerator(cl_ctx, seed=RNG_SEED)
     nsources = lpot_source.density_discr.nnodes
     noise = rng.uniform(queue, nsources, dtype=np.float, a=0.01, b=1.0)
-    panel_sizes = lpot_source.panel_sizes("nsources").with_queue(queue)
+    tunnel_radius = \
+            lpot_source._close_target_tunnel_radius("nsources").with_queue(queue)
 
     def targets_from_sources(sign, dist):
         from pytential import sym, bind
@@ -254,8 +259,8 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
 
     from pytential.target import PointsTarget
 
-    int_targets = PointsTarget(targets_from_sources(-1, noise * panel_sizes / 2))
-    ext_targets = PointsTarget(targets_from_sources(+1, noise * panel_sizes / 2))
+    int_targets = PointsTarget(targets_from_sources(-1, noise * tunnel_radius))
+    ext_targets = PointsTarget(targets_from_sources(+1, noise * tunnel_radius))
     far_targets = PointsTarget(targets_from_sources(+1, FAR_TARGET_DIST_FROM_SOURCE))
 
     # Create target discretizations.
@@ -290,7 +295,7 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
         QBXTargetAssociator(cl_ctx)(lpot_source, target_discrs)
         .get(queue=queue))
 
-    panel_sizes = lpot_source.panel_sizes("nsources").get(queue)
+    expansion_radii = lpot_source._expansion_radii("nsources").get(queue)
 
     int_centers = np.array([axis.get(queue) for axis in int_centers])
     ext_centers = np.array([axis.get(queue) for axis in ext_centers])
@@ -310,7 +315,7 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
                             target_to_source_result, target_to_side_result):
         assert (target_to_side_result == true_side).all()
         dists = la.norm((targets.T - centers.T[target_to_source_result]), axis=1)
-        assert (dists <= panel_sizes[target_to_source_result] / 2).all()
+        assert (dists <= expansion_radii[target_to_source_result]).all()
 
     # Checks that far targets are not assigned a center.
     def check_far_targets(target_to_source_result):
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index 870cc6f0438502e767dc23426d839243619a3023..af43bb94307ebce188706ea072a11ee419f69592 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -160,7 +160,10 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order):
 
         if 0:
             # plot geometry, centers, normals
-            centers = qbx.centers(density_discr, 1)
+
+            from pytential.qbx.utils import get_centers_on_side
+            centers = get_centers_on_side(qbx, 1)
+
             nodes_h = nodes.get()
             centers_h = [centers[0].get(), centers[1].get()]
             pt.plot(nodes_h[0], nodes_h[1], "x-")
diff --git a/test/test_muller.py b/test/test_muller.py
index 30d1a47d8564f607a125bdfd0df5f6e67f811606..fb23e911d32c1c720a3284004014cb83088d75d7 100644
--- a/test/test_muller.py
+++ b/test/test_muller.py
@@ -37,10 +37,14 @@ def test_muller(true_roots):
     :arg n: number of zeros sought
     :return: (roots, niter)
     """
+    np.random.seed(15)
+    z_start = np.random.rand(3) + 1j*np.random.rand(3)
+
     eps = 1e-12
     from pytential.muller import muller_deflate
     roots, niter = muller_deflate(
-        lambda z: poly_with_roots(z, true_roots), len(true_roots))
+        lambda z: poly_with_roots(z, true_roots), len(true_roots),
+        z_start=z_start)
 
     for r_i in roots:
         min_dist, true_root = min(