diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py
index ffc3e0aa9786e1db5b4ec3ab49c86d4e1273cabc..9fd26658a9e91fce4f3c41edf5aa7b698171ac41 100644
--- a/pytential/linalg/proxy.py
+++ b/pytential/linalg/proxy.py
@@ -45,7 +45,6 @@ Proxy Point Generation
 .. autoclass:: ProxyGenerator
 
 .. autofunction:: partition_by_nodes
-.. autofunction:: partition_by_elements
 .. autofunction:: partition_from_coarse
 
 .. autofunction:: gather_block_neighbor_points
@@ -119,85 +118,6 @@ def partition_by_nodes(discr,
                                 ranges.with_queue(None))
 
 
-def partition_by_elements(discr,
-                          use_tree=True,
-                          max_elements_in_box=None):
-    """Generate equally sized ranges of points. The partition is created at the
-    element level, so that all the nodes belonging to an element belong to
-    the same range. This can result in slightly larger differences in size
-    between the ranges, but can be very useful when the individual partitions
-    need to be resampled, integrated, etc.
-
-    :arg discr: a :class:`meshmode.discretization.Discretization`.
-    :arg use_tree: if ``True``, node partitions are generated using a
-        :class:`boxtree.TreeBuilder`, which leads to geometrically close
-        points to belong to the same partition. If ``False``, a simple linear
-        partition is constructed.
-    :arg max_elements_in_box: passed to :class:`boxtree.TreeBuilder`.
-
-    :return: a :class:`sumpy.tools.BlockIndexRanges`.
-    """
-
-    if max_elements_in_box is None:
-        # NOTE: keep in sync with partition_by_nodes
-        max_nodes_in_box = 32
-
-        nunit_nodes = int(np.mean([g.nunit_nodes for g in discr.groups]))
-        max_elements_in_box = max_nodes_in_box // nunit_nodes
-
-    with cl.CommandQueue(discr.cl_context) as queue:
-        if use_tree:
-            from boxtree import box_flags_enum
-            from boxtree import TreeBuilder
-
-            builder = TreeBuilder(discr.cl_context)
-
-            from pytential.qbx.utils import element_centers_of_mass
-            elranges = np.cumsum([group.nelements for group in discr.mesh.groups])
-            elcenters = element_centers_of_mass(discr)
-
-            tree, _ = builder(queue, elcenters,
-                max_particles_in_box=max_elements_in_box)
-
-            groups = discr.groups
-            tree = tree.get(queue)
-            leaf_boxes, = (tree.box_flags
-                           & box_flags_enum.HAS_CHILDREN == 0).nonzero()
-
-            indices = np.empty(len(leaf_boxes), dtype=np.object)
-            for i, ibox in enumerate(leaf_boxes):
-                box_start = tree.box_source_starts[ibox]
-                box_end = box_start + tree.box_source_counts_cumul[ibox]
-
-                ielement = tree.user_source_ids[box_start:box_end]
-                igroup = np.digitize(ielement, elranges)
-
-                indices[i] = np.hstack([_element_node_range(groups[j], k)
-                                        for j, k in zip(igroup, ielement)])
-        else:
-            nelements = discr.mesh.nelements
-            elements = np.array_split(np.arange(0, nelements),
-                                      nelements // max_elements_in_box)
-
-            elranges = np.cumsum([g.nelements for g in discr.groups])
-            elgroups = [np.digitize(elements[i], elranges)
-                        for i in range(len(elements))]
-
-            indices = np.empty(len(elements), dtype=np.object)
-            for i in range(indices.shape[0]):
-                indices[i] = np.hstack([_element_node_range(discr.groups[j], k)
-                                        for j, k in zip(elgroups[i], elements[i])])
-
-        ranges = to_device(queue,
-                np.cumsum([0] + [b.shape[0] for b in indices]))
-        indices = to_device(queue, np.hstack(indices))
-        assert ranges[-1] == discr.nnodes
-
-        return BlockIndexRanges(discr.cl_context,
-                                indices.with_queue(None),
-                                ranges.with_queue(None))
-
-
 def partition_from_coarse(resampler, from_indices):
     """Generate a partition of nodes from an existing partition on a
     coarser discretization. The new partition is generated based on element
@@ -383,7 +303,6 @@ class ProxyGenerator(object):
         radius_expr = radius_expr.format(ratio=self.ratio)
 
         # NOTE: centers of mass are computed using a second-order approximation
-        # that currently matches what is in `element_centers_of_mass`.
         knl = lp.make_kernel([
             "{[irange]: 0 <= irange < nranges}",
             "{[i]: 0 <= i < npoints}",
@@ -410,7 +329,7 @@ class ProxyGenerator(object):
                         (proxy_center[idim, irange] -
                          center_ext[idim, srcindices[i + ioffset]]) ** 2)) + \
                          expansion_radii[srcindices[i + ioffset]])
-                <> rqbx = if(rqbx_ext < rqbx_int, rqbx_int, rqbx_ext)
+                <> rqbx = rqbx_int if rqbx_ext < rqbx_int else rqbx_ext
 
                 proxy_radius[irange] = {radius_expr}
             end
@@ -435,7 +354,6 @@ class ProxyGenerator(object):
             lang_version=MOST_RECENT_LANGUAGE_VERSION)
 
         knl = lp.tag_inames(knl, "idim*:unr")
-
         return knl
 
     @memoize_method
@@ -466,14 +384,20 @@ class ProxyGenerator(object):
         def _affine_map(v, A, b):
             return np.dot(A, v) + b
 
-        from pytential.qbx.utils import get_centers_on_side
+        from pytential import bind, sym
+        radii = bind(self.source,
+                sym.expansion_radii(self.source.ambient_dim))(queue)
+        center_int = bind(self.source,
+                sym.expansion_centers(self.source.ambient_dim, -1))(queue)
+        center_ext = bind(self.source,
+                sym.expansion_centers(self.source.ambient_dim, +1))(queue)
 
         knl = self.get_kernel()
         _, (centers_dev, radii_dev,) = knl(queue,
             sources=self.source.density_discr.nodes(),
-            center_int=get_centers_on_side(self.source, -1),
-            center_ext=get_centers_on_side(self.source, +1),
-            expansion_radii=self.source._expansion_radii("nsources"),
+            center_int=center_int,
+            center_ext=center_ext,
+            expansion_radii=radii,
             srcindices=indices.indices,
             srcranges=indices.ranges, **kwargs)
         centers = centers_dev.get()
diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index bc29d7faec9979b7c9555aab2fd90e68e46947b4..c6762642ce8e4b37c056bd402592fe06b55ee83c 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -56,6 +56,13 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
     .. attribute :: qbx_order
     .. attribute :: fmm_order
 
+    .. automethod :: __init__
+    .. automethod :: with_refinement
+    .. automethod :: copy
+
+    .. attribute :: stage2_density_discr
+    .. attribute :: quad_stage2_density_discr
+
     See :ref:`qbxguts` for some information on the inner workings of this.
     """
 
@@ -348,23 +355,11 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
 
     @memoize_method
     def weights_and_area_elements(self):
-        import pytential.symbolic.primitives as p
-        from pytential.symbolic.execution import bind
+        from pytential import bind, sym
         with cl.CommandQueue(self.cl_context) as queue:
-            # quad_stage2_density_discr is not guaranteed to be usable for
-            # interpolation/differentiation. Use density_discr to find
-            # area element instead, then upsample that.
-
-            area_element = self.refined_interp_to_ovsmp_quad_connection(
-                    queue,
-                    bind(
-                        self.stage2_density_discr,
-                        p.area_element(self.ambient_dim, self.dim)
-                        )(queue))
-
-            qweight = bind(self.quad_stage2_density_discr, p.QWeight())(queue)
-
-            return (area_element.with_queue(queue)*qweight).with_queue(None)
+            return bind(self, sym.weights_and_area_elements(
+                self.ambient_dim,
+                dofdesc=sym.QBX_SOURCE_QUAD_STAGE2))(queue).with_queue(None)
 
     # }}}
 
@@ -389,7 +384,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
         .. warning::
 
             This always returns a
-            :class:`~meshmode.discretization.connection.DirectDiscretizationConnect`.
+            :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`.
             In case the geometry has been refined multiple times, a direct
             connection can have a large number of groups and/or
             interpolation batches, making it scale significantly worse than
@@ -465,114 +460,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
 
         return lpot, connection
 
-    @property
-    @memoize_method
-    def h_max(self):
-        with cl.CommandQueue(self.cl_context) as queue:
-            quad_res = self._coarsest_quad_resolution("npanels").with_queue(queue)
-            return cl.array.max(quad_res).get().item()
-
     # {{{ internal API
 
-    @memoize_method
-    def _panel_centers_of_mass(self):
-        import pytential.qbx.utils as utils
-        return utils.element_centers_of_mass(self.density_discr)
-
-    @memoize_method
-    def _stage2_panel_centers_of_mass(self):
-        import pytential.qbx.utils as utils
-        return utils.element_centers_of_mass(self.stage2_density_discr)
-
-    def _dim_fudge_factor(self):
-        if self.density_discr.dim == 2:
-            return 0.5
-        else:
-            return 1
-
-    @memoize_method
-    def _expansion_radii(self, last_dim_length):
-        with cl.CommandQueue(self.cl_context) as queue:
-            return (self._coarsest_quad_resolution(last_dim_length)
-                    .with_queue(queue)
-                    * 0.5 * self._dim_fudge_factor()).with_queue(None)
-
-    # _expansion_radii should not be needed for the fine discretization
-
-    @memoize_method
-    def _source_danger_zone_radii(self, last_dim_length="npanels"):
-        # This should be the expression of the expansion radii, but
-        #
-        # - in reference to the stage 2 discretization
-        # - mutliplied by 0.75 because
-        #
-        #   - Setting this equal to the expansion radii ensures that *every*
-        #     stage 2 element will be refined, which is wasteful.
-        #     (so this needs to be smaller than that)
-        #
-
-        #   - Setting this equal to half the expansion radius will not provide
-        #     a refinement 'buffer layer' at a 2x coarsening fringe.
-
-        with cl.CommandQueue(self.cl_context) as queue:
-            return (
-                    (self._stage2_coarsest_quad_resolution(last_dim_length)
-                        .with_queue(queue))
-                    * 0.5 * 0.75 * self._dim_fudge_factor()).with_queue(None)
-
-    @memoize_method
-    def _close_target_tunnel_radius(self, last_dim_length):
-        with cl.CommandQueue(self.cl_context) as queue:
-            return (
-                    self._expansion_radii(last_dim_length).with_queue(queue)
-                    * 0.5
-                    ).with_queue(None)
-
-    @memoize_method
-    def _coarsest_quad_resolution(self, last_dim_length="npanels"):
-        """This measures the quadrature resolution across the
-        mesh. In a 1D uniform mesh of uniform 'parametrization speed', it
-        should be the same as the panel length.
-        """
-        import pytential.qbx.utils as utils
-        from pytential import sym, bind
-        with cl.CommandQueue(self.cl_context) as queue:
-            maxstretch = bind(
-                    self,
-                    sym._simplex_mapping_max_stretch_factor(
-                        self.ambient_dim)
-                    )(queue)
-
-            maxstretch = utils.to_last_dim_length(
-                    self.density_discr, maxstretch, last_dim_length)
-            maxstretch = maxstretch.with_queue(None)
-
-        return maxstretch
-
-    @memoize_method
-    def _stage2_coarsest_quad_resolution(self, last_dim_length="npanels"):
-        """This measures the quadrature resolution across the
-        mesh. In a 1D uniform mesh of uniform 'parametrization speed', it
-        should be the same as the panel length.
-        """
-        if last_dim_length != "npanels":
-            # Not technically required below, but no need to loosen for now.
-            raise NotImplementedError()
-
-        import pytential.qbx.utils as utils
-        from pytential import sym, bind
-        with cl.CommandQueue(self.cl_context) as queue:
-            maxstretch = bind(
-                    self, sym._simplex_mapping_max_stretch_factor(
-                        self.ambient_dim,
-                        where=sym.QBXSourceStage2(sym.DEFAULT_SOURCE))
-                    )(queue)
-            maxstretch = utils.to_last_dim_length(
-                    self.stage2_density_discr, maxstretch, last_dim_length)
-            maxstretch = maxstretch.with_queue(None)
-
-        return maxstretch
-
     @memoize_method
     def qbx_fmm_geometry_data(self, target_discrs_and_qbx_sides):
         """
@@ -687,7 +576,7 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
                 tgt_name_and_side_to_number[key] = \
                         len(target_discrs_and_qbx_sides)
 
-                target_discr = bound_expr.places[o.target_name]
+                target_discr = bound_expr.places.get_geometry(o.target_name)
                 if isinstance(target_discr, LayerPotentialSourceBase):
                     target_discr = target_discr.density_discr
 
@@ -838,7 +727,15 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
         strengths = (evaluate(insn.density).with_queue(queue)
                 * self.weights_and_area_elements())
 
-        import pytential.qbx.utils as utils
+        from pytential import bind, sym
+        expansion_radii = bind(self,
+                sym.expansion_radii(self.ambient_dim))(queue)
+        centers = {
+                -1: bind(self,
+                    sym.expansion_centers(self.ambient_dim, -1))(queue),
+                +1: bind(self,
+                    sym.expansion_centers(self.ambient_dim, +1))(queue)
+                }
 
         # FIXME: Do this all at once
         result = []
@@ -855,9 +752,9 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
                 evt, output_for_each_kernel = lpot_applier(
                         queue, target_discr.nodes(),
                         self.quad_stage2_density_discr.nodes(),
-                        utils.get_centers_on_side(self, o.qbx_forced_limit),
+                        centers[o.qbx_forced_limit],
                         [strengths],
-                        expansion_radii=self._expansion_radii("nsources"),
+                        expansion_radii=expansion_radii,
                         **kernel_args)
                 result.append((o.name, output_for_each_kernel[o.kernel_index]))
             else:
diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py
index a80a8098aa234834c216aff5cdf0e5b3e7ccd39e..829e706c481227bbc49b3e524094555cb1f541b8 100644
--- a/pytential/qbx/refinement.py
+++ b/pytential/qbx/refinement.py
@@ -304,7 +304,10 @@ class RefinerWrangler(TreeWranglerBase):
         found_panel_to_refine.finish()
         unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
 
-        center_danger_zone_radii = lpot_source._expansion_radii("ncenters")
+        from pytential import bind, sym
+        center_danger_zone_radii = bind(lpot_source, sym.expansion_radii(
+            lpot_source.ambient_dim,
+            granularity=sym.GRANULARITY_CENTER))(self.queue)
 
         evt = knl(
             *unwrap_args(
@@ -358,8 +361,11 @@ class RefinerWrangler(TreeWranglerBase):
         found_panel_to_refine = cl.array.zeros(self.queue, 1, np.int32)
         found_panel_to_refine.finish()
 
-        source_danger_zone_radii_by_panel = \
-                lpot_source._source_danger_zone_radii("npanels")
+        from pytential import bind, sym
+        source_danger_zone_radii_by_panel = bind(lpot_source,
+                sym._source_danger_zone_radii(
+                    lpot_source.ambient_dim,
+                    dofdesc=sym.GRANULARITY_ELEMENT))(self.queue)
         unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
 
         evt = knl(
@@ -399,7 +405,6 @@ class RefinerWrangler(TreeWranglerBase):
 
         evt, out = knl(self.queue,
                        element_property=element_property,
-                       # lpot_source._coarsest_quad_resolution("npanels")),
                        refine_flags=refine_flags,
                        refine_flags_updated=np.array(0),
                        threshold=np.array(threshold),
@@ -599,11 +604,14 @@ def refine_for_global_qbx(lpot_source, wrangler,
             with ProcessLogger(logger,
                     "checking kernel length scale to panel size ratio"):
 
+                from pytential import bind, sym
+                quad_resolution = bind(lpot_source, sym._quad_resolution(
+                    lpot_source.ambient_dim,
+                    dofdesc=sym.GRANULARITY_ELEMENT))(wrangler.queue)
+
                 violates_kernel_length_scale = \
                         wrangler.check_element_prop_threshold(
-                                element_property=(
-                                    lpot_source._coarsest_quad_resolution(
-                                        "npanels")),
+                                element_property=quad_resolution,
                                 threshold=kernel_length_scale,
                                 refine_flags=refine_flags, debug=debug)
 
@@ -615,15 +623,11 @@ def refine_for_global_qbx(lpot_source, wrangler,
         if scaled_max_curvature_threshold is not None:
             with ProcessLogger(logger,
                     "checking scaled max curvature threshold"):
-                from pytential.qbx.utils import to_last_dim_length
                 from pytential import sym, bind
-                scaled_max_curv = to_last_dim_length(
-                        lpot_source.density_discr,
-                        bind(lpot_source,
-                            sym.ElementwiseMax(
-                                sym._scaled_max_curvature(
-                                    lpot_source.density_discr.ambient_dim)))
-                            (wrangler.queue), "npanels")
+                scaled_max_curv = bind(lpot_source,
+                    sym.ElementwiseMax(
+                        sym._scaled_max_curvature(lpot_source.ambient_dim),
+                        dofdesc=sym.GRANULARITY_ELEMENT))(wrangler.queue)
 
                 violates_scaled_max_curv = \
                         wrangler.check_element_prop_threshold(
diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py
index 25170c5f593f136576005937dfb4f838264606a4..01fd6eb653cbeda02d812589638f98ee4b0abfe5 100644
--- a/pytential/qbx/target_assoc.py
+++ b/pytential/qbx/target_assoc.py
@@ -452,12 +452,12 @@ class TargetAssociationWrangler(TreeWranglerBase):
         found_target_close_to_panel.finish()
 
         # Perform a space invader query over the sources.
+        from pytential import bind, sym
         source_slice = tree.sorted_target_ids[tree.qbx_user_source_slice]
         sources = [
                 axis.with_queue(self.queue)[source_slice] for axis in tree.sources]
-        tunnel_radius_by_source = (
-                lpot_source._close_target_tunnel_radius("nsources")
-                .with_queue(self.queue))
+        tunnel_radius_by_source = bind(lpot_source,
+                sym._close_target_tunnel_radii(lpot_source.ambient_dim))(self.queue)
 
         # Target-marking algorithm (TGTMARK):
         #
@@ -493,7 +493,8 @@ class TargetAssociationWrangler(TreeWranglerBase):
                 wait_for=wait_for)
         wait_for = [evt]
 
-        tunnel_radius_by_source = lpot_source._close_target_tunnel_radius("nsources")
+        tunnel_radius_by_source = bind(lpot_source,
+            sym._close_target_tunnel_radii(lpot_source.ambient_dim))(self.queue)
 
         evt = knl(
             *unwrap_args(
@@ -544,13 +545,15 @@ class TargetAssociationWrangler(TreeWranglerBase):
             marked_target_count = int(cl.array.sum(target_status).get())
 
         # Perform a space invader query over the centers.
+        from pytential import bind, sym
         center_slice = (
                 tree.sorted_target_ids[tree.qbx_user_center_slice]
                 .with_queue(self.queue))
         centers = [
                 axis.with_queue(self.queue)[center_slice] for axis in tree.sources]
-        expansion_radii_by_center = \
-                lpot_source._expansion_radii("ncenters").with_queue(self.queue)
+        expansion_radii_by_center = bind(lpot_source, sym.expansion_radii(
+            lpot_source.ambient_dim,
+            granularity=sym.GRANULARITY_CENTER))(self.queue)
         expansion_radii_by_center_with_tolerance = \
                 expansion_radii_by_center * (1 + target_association_tolerance)
 
@@ -625,12 +628,12 @@ class TargetAssociationWrangler(TreeWranglerBase):
         found_panel_to_refine.finish()
 
         # Perform a space invader query over the sources.
+        from pytential import bind, sym
         source_slice = tree.user_source_ids[tree.qbx_user_source_slice]
         sources = [
                 axis.with_queue(self.queue)[source_slice] for axis in tree.sources]
-        tunnel_radius_by_source = (
-                lpot_source._close_target_tunnel_radius("nsources")
-                .with_queue(self.queue))
+        tunnel_radius_by_source = bind(lpot_source,
+                sym._close_target_tunnel_radii(lpot_source.ambient_dim))(self.queue)
 
         # See (TGTMARK) above for algorithm.
 
@@ -643,6 +646,9 @@ class TargetAssociationWrangler(TreeWranglerBase):
                 wait_for=wait_for)
         wait_for = [evt]
 
+        tunnel_radius_by_source = bind(lpot_source,
+                sym._close_target_tunnel_radii(lpot_source.ambient_dim))(self.queue)
+
         evt = knl(
             *unwrap_args(
                 tree, peer_lists,
@@ -653,7 +659,7 @@ class TargetAssociationWrangler(TreeWranglerBase):
                 tree.qbx_user_target_slice.start,
                 tree.nqbxpanels,
                 tree.sorted_target_ids,
-                lpot_source._close_target_tunnel_radius("nsources"),
+                tunnel_radius_by_source,
                 target_status,
                 box_to_search_dist,
                 refine_flags,
diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py
index 01554e0f322f8a4805c1bec79198ee0d801904b9..194c90fc7eb9ff18a2d3c6bfb1da6df33a3f1951 100644
--- a/pytential/qbx/utils.py
+++ b/pytential/qbx/utils.py
@@ -26,13 +26,11 @@ THE SOFTWARE.
 """
 
 
-import loopy as lp
 import numpy as np
 from boxtree.tree import Tree
 import pyopencl as cl
 import pyopencl.array # noqa
-from pytools import memoize, memoize_method
-from loopy.version import MOST_RECENT_LANGUAGE_VERSION
+from pytools import memoize_method
 
 import logging
 logger = logging.getLogger(__name__)
@@ -47,7 +45,6 @@ QBX_TREE_C_PREAMBLE = r"""//CL:mako//
 // have their own numbering starting at 0. These macros convert
 // the per-class numbering into the internal tree particle number.
 #define INDEX_FOR_CENTER_PARTICLE(i) (sorted_target_ids[center_offset + i])
-#define INDEX_FOR_PANEL_PARTICLE(i) (sorted_target_ids[panel_offset + i])
 #define INDEX_FOR_SOURCE_PARTICLE(i) (sorted_target_ids[source_offset + i])
 #define INDEX_FOR_TARGET_PARTICLE(i) (sorted_target_ids[target_offset + i])
 
@@ -70,31 +67,6 @@ QBX_TREE_MAKO_DEFS = r"""//CL:mako//
 # }}}
 
 
-# {{{ interleaver kernel
-
-@memoize
-def get_interleaver_kernel(dtype):
-    # NOTE: Returned kernel needs dstlen or dst parameter
-    from pymbolic import var
-    knl = lp.make_kernel(
-        "[srclen,dstlen] -> {[i]: 0<=i<srclen}",
-        """
-        dst[2*i] = src1[i]
-        dst[2*i+1] = src2[i]
-        """, [
-            lp.GlobalArg("src1", shape=(var("srclen"),), dtype=dtype),
-            lp.GlobalArg("src2", shape=(var("srclen"),), dtype=dtype),
-            lp.GlobalArg("dst", shape=(var("dstlen"),), dtype=dtype),
-            "..."
-        ],
-        assumptions="2*srclen = dstlen",
-        lang_version=MOST_RECENT_LANGUAGE_VERSION)
-    knl = lp.split_iname(knl, "i", 128, inner_tag="l.0", outer_tag="g.0")
-    return knl
-
-# }}}
-
-
 # {{{ make interleaved centers
 
 def get_interleaved_centers(queue, lpot_source):
@@ -102,22 +74,15 @@ def get_interleaved_centers(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)
-    int_centers = get_centers_on_side(lpot_source, -1)
-    ext_centers = get_centers_on_side(lpot_source, +1)
-
-    result = []
-    wait_for = []
-
-    for int_axis, ext_axis in zip(int_centers, ext_centers):
-        axis = cl.array.empty(queue, len(int_axis) * 2, int_axis.dtype)
-        evt, _ = knl(queue, src1=int_axis, src2=ext_axis, dst=axis)
-        result.append(axis)
-        wait_for.append(evt)
+    from pytential import bind, sym
+    int_centers = bind(lpot_source,
+            sym.expansion_centers(lpot_source.ambient_dim, -1))(queue)
+    ext_centers = bind(lpot_source,
+            sym.expansion_centers(lpot_source.ambient_dim, +1))(queue)
 
-    cl.wait_for_events(wait_for)
-
-    return result
+    from pytential.symbolic.dof_connection import CenterGranularityConnection
+    interleaver = CenterGranularityConnection(lpot_source.density_discr)
+    return interleaver(queue, [int_centers, ext_centers])
 
 # }}}
 
@@ -129,14 +94,14 @@ 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")
+    from pytential import bind, sym
 
-    result = cl.array.empty(queue, len(radii) * 2, radii.dtype)
-    evt, _ = knl(queue, src1=radii, src2=radii, dst=result)
-    evt.wait()
+    radii = bind(lpot_source,
+            sym.expansion_radii(lpot_source.ambient_dim))(queue)
 
-    return result
+    from pytential.symbolic.dof_connection import CenterGranularityConnection
+    interleaver = CenterGranularityConnection(lpot_source.density_discr)
+    return interleaver(queue, radii)
 
 # }}}
 
@@ -211,126 +176,6 @@ class TreeWranglerBase(object):
 # }}}
 
 
-# {{{ to_last_dim_length
-
-def to_last_dim_length(discr, vec, last_dim_length, queue=None):
-    """Takes a :class:`pyopencl.array.Array` with a last axis that has the same
-    length as the number of discretization nodes in the discretization *discr*
-    and converts it so that the last axis has a length as specified by
-    *last_dim_length*.
-    """
-
-    queue = queue or vec.queue
-
-    if last_dim_length == "nsources":
-        return vec
-
-    elif last_dim_length == "ncenters":
-        knl = get_interleaver_kernel(vec.dtype)
-        _, (result,) = knl(queue, dstlen=2*discr.nnodes, src1=vec, src2=vec)
-        return result
-
-    elif last_dim_length == "npanels":
-        knl = lp.make_kernel(
-            "{[i,k]: 0<=i<nelements}",
-            "result[i] = a[i,0]",
-            [
-                lp.GlobalArg("a", shape="nelements, nunit_nodes", dtype=None),
-                lp.ValueArg("nunit_nodes", dtype=np.int32),
-                "..."
-                ],
-            name="subsample_to_panels",
-            lang_version=MOST_RECENT_LANGUAGE_VERSION,
-            )
-
-        knl = lp.split_iname(knl, "i", 128, inner_tag="l.0", outer_tag="g.0")
-
-        result = cl.array.empty(queue, discr.mesh.nelements, vec.dtype)
-        for group_nr, group in enumerate(discr.groups):
-            knl(queue,
-                nelements=group.nelements,
-                a=group.view(vec),
-                result=el_view(discr, group_nr, result))
-
-        return result
-    else:
-        raise ValueError("unknown dim length specified")
-
-# }}}
-
-
-# {{{ 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",
-        lang_version=MOST_RECENT_LANGUAGE_VERSION)
-
-    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.mesh.groups[group_nr]
-
-    return global_array[
-        ..., group.element_nr_base:group.element_nr_base + group.nelements] \
-        .reshape(
-            global_array.shape[:-1]
-            + (group.nelements,))
-
-# }}}
-
-
 # {{{ tree-with-metadata: data structure
 
 class TreeWithQBXMetadata(Tree):
@@ -404,7 +249,6 @@ class TreeWithQBXMetadata(Tree):
 
     .. attribute:: qbx_user_source_slice
     .. attribute:: qbx_user_center_slice
-    .. attribute:: qbx_user_panel_slice
     .. attribute:: qbx_user_target_slice
     """
     pass
@@ -424,7 +268,7 @@ def build_tree_with_qbx_metadata(
     """Return a :class:`TreeWithQBXMetadata` built from the given layer
     potential source. This contains particles of four different types:
 
-       * source particles and panel centers of mass either from
+       * source particles either from
          ``lpot_source.density_discr`` or
          ``lpot_source.stage2_density_discr``
        * centers from ``lpot_source.centers()``
@@ -437,37 +281,31 @@ def build_tree_with_qbx_metadata(
 
     :arg targets_list: A list of :class:`pytential.target.TargetBase`
 
-    :arg use_stage2_discr: If *True*, builds a tree with sources/centers of
-        mass from ``lpot_source.stage2_density_discr``. If *False* (default),
+    :arg use_stage2_discr: If *True*, builds a tree with sources
+        from ``lpot_source.stage2_density_discr``. If *False* (default),
         they are from ``lpot_source.density_discr``.
     """
     # The ordering of particles is as follows:
     # - sources go first
     # - then centers
-    # - then panels (=centers of mass)
     # - then targets
 
-    sources = (
-            lpot_source.density_discr.nodes()
-            if not use_stage2_discr
-            else lpot_source.quad_stage2_density_discr.nodes())
+    if use_stage2_discr:
+        density_discr = lpot_source.quad_stage2_density_discr
+    else:
+        density_discr = lpot_source.density_discr
 
+    sources = density_discr.nodes()
     centers = get_interleaved_centers(queue, lpot_source)
-
-    centers_of_mass = (
-            lpot_source._panel_centers_of_mass()
-            if not use_stage2_discr
-            else lpot_source._stage2_panel_centers_of_mass())
-
     targets = (tgt.nodes() for tgt in targets_list)
 
     particles = tuple(
             cl.array.concatenate(dim_coords, queue=queue)
-            for dim_coords in zip(sources, centers, centers_of_mass, *targets))
+            for dim_coords in zip(sources, centers, *targets))
 
     # Counts
     nparticles = len(particles[0])
-    npanels = len(centers_of_mass[0])
+    npanels = density_discr.mesh.nelements
     nsources = len(sources[0])
     ncenters = len(centers[0])
     # Each source gets an interior / exterior center.
@@ -481,12 +319,10 @@ def build_tree_with_qbx_metadata(
     qbx_user_center_slice = slice(center_slice_start, center_slice_start + ncenters)
 
     panel_slice_start = center_slice_start + ncenters
-    qbx_user_panel_slice = slice(panel_slice_start, panel_slice_start + npanels)
-
-    target_slice_start = panel_slice_start + npanels
+    target_slice_start = panel_slice_start
     qbx_user_target_slice = slice(target_slice_start, target_slice_start + ntargets)
 
-    # Build tree with sources, centers, and centers of mass. Split boxes
+    # Build tree with sources and centers. Split boxes
     # only because of sources.
     refine_weights = cl.array.zeros(queue, nparticles, np.int32)
     refine_weights[:nsources].fill(1)
@@ -505,8 +341,7 @@ def build_tree_with_qbx_metadata(
     for class_name, particle_slice, fixup in (
             ("box_to_qbx_source", qbx_user_source_slice, 0),
             ("box_to_qbx_target", qbx_user_target_slice, -target_slice_start),
-            ("box_to_qbx_center", qbx_user_center_slice, -center_slice_start),
-            ("box_to_qbx_panel", qbx_user_panel_slice, -panel_slice_start)):
+            ("box_to_qbx_center", qbx_user_center_slice, -center_slice_start)):
         flags.fill(0)
         flags[particle_slice].fill(1)
         flags.finish()
@@ -525,11 +360,6 @@ def build_tree_with_qbx_metadata(
     del box_to_class
 
     # Compute panel => source relation
-    if use_stage2_discr:
-        density_discr = lpot_source.quad_stage2_density_discr
-    else:
-        density_discr = lpot_source.density_discr
-
     qbx_panel_to_source_starts = cl.array.empty(
             queue, npanels + 1, dtype=tree.particle_id_dtype)
     el_offset = 0
@@ -562,7 +392,6 @@ def build_tree_with_qbx_metadata(
         qbx_panel_to_source_starts=qbx_panel_to_source_starts,
         qbx_panel_to_center_starts=qbx_panel_to_center_starts,
         qbx_user_source_slice=qbx_user_source_slice,
-        qbx_user_panel_slice=qbx_user_panel_slice,
         qbx_user_center_slice=qbx_user_center_slice,
         qbx_user_target_slice=qbx_user_target_slice,
         nqbxpanels=npanels,
diff --git a/pytential/source.py b/pytential/source.py
index e628fe21291341b15e8d24b9a64c858d2b94e754..1a55878ef41a4d9b507c2d54562428bd723ea2c4 100644
--- a/pytential/source.py
+++ b/pytential/source.py
@@ -130,8 +130,7 @@ class PointPotentialSource(PotentialSource):
         for arg_name, arg_expr in six.iteritems(insn.kernel_arguments):
             kernel_args[arg_name] = evaluate(arg_expr)
 
-        strengths = (evaluate(insn.density).with_queue(queue)
-                * self.weights_and_area_elements())
+        strengths = evaluate(insn.density).with_queue(queue).copy()
 
         # FIXME: Do this all at once
         result = []
@@ -153,8 +152,7 @@ class PointPotentialSource(PotentialSource):
     @memoize_method
     def weights_and_area_elements(self):
         with cl.CommandQueue(self.cl_context) as queue:
-            result = cl.array.empty(queue, self._nodes.shape[-1],
-                    dtype=self.real_dtype)
+            result = cl.array.empty(queue, self.nnodes, dtype=self.real_dtype)
             result.fill(1)
 
         return result.with_queue(None)
@@ -174,7 +172,6 @@ class LayerPotentialSourceBase(PotentialSource):
     .. attribute:: stage2_density_discr
     .. attribute:: quad_stage2_density_discr
     .. attribute:: resampler
-    .. method:: with_refinement
 
     .. rubric:: Discretization data
 
@@ -183,7 +180,6 @@ class LayerPotentialSourceBase(PotentialSource):
     .. attribute:: dim
     .. attribute:: real_dtype
     .. attribute:: complex_dtype
-    .. attribute:: h_max
 
     .. rubric:: Execution
 
diff --git a/pytential/symbolic/compiler.py b/pytential/symbolic/compiler.py
index cb9a9c1de493da6df5b005c8e8d665d608e40da4..8f7ca490a22ccccc3c4e45a185ab9c596d429d31 100644
--- a/pytential/symbolic/compiler.py
+++ b/pytential/symbolic/compiler.py
@@ -224,7 +224,7 @@ class ComputePotentialInstruction(Instruction):
                 ", ".join(args), "\n  ".join(lines))
 
     def get_exec_function(self, exec_mapper):
-        source = exec_mapper.bound_expr.places[self.source]
+        source = exec_mapper.bound_expr.places.get_geometry(self.source)
         return source.exec_compute_potential_insn
 
     def __hash__(self):
@@ -440,8 +440,9 @@ class OperatorCompiler(IdentityMapper):
 
     def op_group_features(self, expr):
         from pytential.symbolic.primitives import hashable_kernel_args
+        lpot_source = self.places.get_geometry(expr.source)
         return (
-                self.places[expr.source].op_group_features(expr)
+                lpot_source.op_group_features(expr)
                 + hashable_kernel_args(expr.kernel_arguments))
 
     @memoize_method
diff --git a/pytential/symbolic/dof_connection.py b/pytential/symbolic/dof_connection.py
new file mode 100644
index 0000000000000000000000000000000000000000..8cc5d1e4e8202dee9a20fed4013f4f106fd75a60
--- /dev/null
+++ b/pytential/symbolic/dof_connection.py
@@ -0,0 +1,266 @@
+# -*- coding: utf-8 -*-
+from __future__ import division, absolute_import, print_function
+
+__copyright__ = """
+Copyright (C) 2016 Matt Wala
+Copyright (C) 2019 Alexandru Fikl
+"""
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import pyopencl as cl
+import pyopencl.array # noqa
+from pytools import memoize
+
+import loopy as lp
+from loopy.version import MOST_RECENT_LANGUAGE_VERSION
+
+
+__doc__ = """
+
+Connections
+-----------
+
+.. autoclass:: GranularityConnection
+.. autoclass:: CenterGranularityConnection
+.. autoclass:: DOFConnection
+.. autofunction:: connection_from_dds
+
+"""
+
+
+# {{{ granularity connections
+
+class GranularityConnection(object):
+    """Abstract interface for transporting a DOF between different levels
+    of granularity.
+
+    .. attribute:: discr
+    .. automethod:: __call__
+    """
+
+    def __init__(self, discr):
+        self.discr = discr
+
+    @property
+    def from_discr(self):
+        return self.discr
+
+    @property
+    def to_discr(self):
+        return self.discr
+
+    def __call__(self, queue, vec):
+        raise NotImplementedError()
+
+
+class CenterGranularityConnection(GranularityConnection):
+    """A :class:`GranularityConnection` used to transport from node data
+    (:class:`~pytential.symbolic.primitives.GRANULARITY_NODE`) to expansion
+    centers (:class:`~pytential.symbolic.primitives.GRANULARITY_CENTER`).
+
+    .. attribute:: discr
+    .. automethod:: __call__
+    """
+
+    def __init__(self, discr):
+        super(CenterGranularityConnection, self).__init__(discr)
+
+    @memoize
+    def kernel(self):
+        knl = lp.make_kernel(
+            "[srclen, dstlen] -> {[i]: 0 <= i < srclen}",
+            """
+            dst[2*i] = src1[i]
+            dst[2*i + 1] = src2[i]
+            """,
+            [
+                lp.GlobalArg("src1", shape="srclen"),
+                lp.GlobalArg("src2", shape="srclen"),
+                lp.GlobalArg("dst", shape="dstlen"),
+                "..."
+            ],
+            name="node_interleaver_knl",
+            assumptions="2*srclen = dstlen",
+            lang_version=MOST_RECENT_LANGUAGE_VERSION,
+            )
+
+        knl = lp.split_iname(knl, "i", 128,
+                inner_tag="l.0", outer_tag="g.0")
+        return knl
+
+    def __call__(self, queue, vecs):
+        r"""
+        :arg vecs: a single :class:`pyopencl.array.Array` or a pair of arrays.
+        :return: an interleaved array or list of :class:`pyopencl.array.Array`s.
+            If *vecs* was a pair of arrays :math:`(x, y)`, they are
+            interleaved as :math:`[x_1, y_1, x_2, y_2, \ddots, x_n, y_n]`.
+            A single array is simply interleaved with itself.
+        """
+
+        if isinstance(vecs, cl.array.Array):
+            vecs = [[vecs], [vecs]]
+        elif isinstance(vecs, (list, tuple)):
+            assert len(vecs) == 2
+        else:
+            raise ValueError('cannot interleave arrays')
+
+        result = []
+        for src1, src2 in zip(vecs[0], vecs[1]):
+            if not isinstance(src1, cl.array.Array) \
+                    or not isinstance(src2, cl.array.Array):
+                raise TypeError('non-array passed to connection')
+
+            if src1.shape != (self.discr.nnodes,) \
+                    or src2.shape != (self.discr.nnodes,):
+                raise ValueError('invalid shape of incoming array')
+
+            axis = cl.array.empty(queue, 2 * len(src1), src1.dtype)
+            self.kernel()(queue,
+                    src1=src1, src2=src2, dst=axis)
+            result.append(axis)
+
+        return result[0] if len(result) == 1 else result
+
+# }}}
+
+
+# {{{ dof connection
+
+class DOFConnection(object):
+    """An interpolation operation for converting a DOF vector between
+    different DOF types, as described by
+    :class:`~pytential.symbolic.primitives.DOFDescriptor`.
+
+    .. attribute:: connections
+
+        A list of
+        :class:`~meshmode.discretization.connection.DiscretizationConnection`s
+        and :class:`GranularityConnection`s used to transport from the given
+        source to the target.
+
+    .. attribute:: from_dd
+
+        A :class:`~pytential.symbolic.primitives.DOFDescriptor` for the
+        DOF type of the incoming array.
+
+    .. attribute:: to_dd
+
+        A :class:`~pytential.symbolic.primitives.DOFDescriptor` for the
+        DOF type of the outgoing array.
+
+    .. attribute:: from_discr
+    .. attribute:: to_discr
+
+    .. automethod:: __call__
+    """
+
+    def __init__(self, connections, from_dd=None, to_dd=None):
+        self.from_dd = from_dd
+        self.to_dd = to_dd
+        self.connections = connections
+
+        from meshmode.discretization.connection import DiscretizationConnection
+        for conn in self.connections:
+            if not isinstance(conn,
+                    (DiscretizationConnection, GranularityConnection)):
+                raise ValueError('unsupported connection type: {}'
+                        .format(type(conn)))
+
+        if self.connections:
+            self.from_discr = self.connections[0].from_discr
+            self.to_discr = self.connections[-1].to_discr
+
+    def __call__(self, queue, vec):
+        for conn in self.connections:
+            vec = conn(queue, vec)
+
+        return vec
+
+
+def connection_from_dds(places, from_dd, to_dd):
+    """
+    :arg places: a :class:`~pytential.symbolic.execution.GeometryCollection`
+        or an argument taken by its constructor.
+    :arg from_dd: a descriptor for the incoming degrees of freedom. This
+        can be a :class:`~pytential.symbolic.primitives.DOFDescriptor`
+        or an identifier that can be transformed into one by
+        :func:`~pytential.symbolic.primitives.as_dofdesc`.
+    :arg to_dd: a descriptor for the outgoing degrees of freedom.
+
+    :return: a :class:`DOFConnection` transporting between the two
+        kinds of DOF vectors.
+    """
+
+    from pytential import sym
+    from_dd = sym.as_dofdesc(from_dd)
+    to_dd = sym.as_dofdesc(to_dd)
+
+    from pytential.symbolic.execution import GeometryCollection
+    if not isinstance(places, GeometryCollection):
+        places = GeometryCollection(places)
+    from_discr = places.get_geometry(from_dd)
+
+    if from_dd.geometry != to_dd.geometry:
+        raise ValueError("cannot interpolate between different geometries")
+
+    if from_dd.granularity is not sym.GRANULARITY_NODE:
+        raise ValueError("can only interpolate from `GRANULARITY_NODE`")
+
+    connections = []
+    if from_dd.discr_stage is not to_dd.discr_stage:
+        from pytential.qbx import QBXLayerPotentialSource
+        if not isinstance(from_discr, QBXLayerPotentialSource):
+            raise ValueError("can only interpolate on a "
+                    "`QBXLayerPotentialSource`")
+
+        if to_dd.discr_stage is not sym.QBX_SOURCE_QUAD_STAGE2:
+            # TODO: can probably extend this to project from a QUAD_STAGE2
+            # using L2ProjectionInverseDiscretizationConnection
+            raise ValueError("can only interpolate to "
+                "`QBX_SOURCE_QUAD_STAGE2`")
+
+        if from_dd.discr_stage is sym.QBX_SOURCE_QUAD_STAGE2:
+            pass
+        elif from_dd.discr_stage is sym.QBX_SOURCE_STAGE2:
+            connections.append(
+                    from_discr.refined_interp_to_ovsmp_quad_connection)
+        else:
+            connections.append(from_discr.resampler)
+
+    if from_dd.granularity is not to_dd.granularity:
+        to_discr = places.get_discretization(to_dd)
+
+        if to_dd.granularity is sym.GRANULARITY_NODE:
+            pass
+        elif to_dd.granularity is sym.GRANULARITY_CENTER:
+            connections.append(CenterGranularityConnection(to_discr))
+        elif to_dd.granularity is sym.GRANULARITY_ELEMENT:
+            raise ValueError("Creating a connection to element granularity "
+                    "is not allowed. Use Elementwise{Max,Min,Sum}.")
+        else:
+            raise ValueError("invalid to_dd granularity: %s" % to_dd.granularity)
+
+    return DOFConnection(connections, from_dd=from_dd, to_dd=to_dd)
+
+# }}}
+
+# vim: foldmethod=marker
diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py
index 34697a6de6ad2a8a00a3782205e8793032351b03..5a8f8b786bd41209700a46a4cefbfee3cdf9a0ad 100644
--- a/pytential/symbolic/execution.py
+++ b/pytential/symbolic/execution.py
@@ -46,6 +46,22 @@ from pytential import sym
 
 # {{{ evaluation mapper
 
+def mesh_el_view(mesh, group_nr, global_array):
+    """Return a view of *global_array* of shape
+    ``(..., mesh.groups[group_nr].nelements)``
+    where *global_array* is of shape ``(..., nelements)``,
+    where *nelements* is the global (per-mesh) element count.
+    """
+
+    group = mesh.groups[group_nr]
+
+    return global_array[
+        ..., group.element_nr_base:group.element_nr_base + group.nelements] \
+        .reshape(
+            global_array.shape[:-1]
+            + (group.nelements,))
+
+
 class EvaluationMapper(EvaluationMapperBase):
     def __init__(self, bound_expr, queue, context={},
             target_geometry=None,
@@ -53,6 +69,7 @@ class EvaluationMapper(EvaluationMapperBase):
         EvaluationMapperBase.__init__(self, context)
 
         self.bound_expr = bound_expr
+        self.places = bound_expr.places
         self.queue = queue
 
     # {{{ map_XXX
@@ -84,30 +101,63 @@ class EvaluationMapper(EvaluationMapperBase):
         return cl.array.max(self.rec(expr.operand)).get()[()]
 
     def _map_elementwise_reduction(self, reduction_name, expr):
-        @memoize_in(self.bound_expr, "elementwise_"+reduction_name)
-        def knl():
+        @memoize_in(self.places, "elementwise_node_"+reduction_name)
+        def node_knl():
             import loopy as lp
             knl = lp.make_kernel(
-                "{[el, idof, jdof]: 0<=el<nelements and 0<=idof,jdof<ndofs}",
-                "result[el, idof] = %s(jdof, input[el, jdof])" % reduction_name,
-                default_offset=lp.auto,
-                lang_version=MOST_RECENT_LANGUAGE_VERSION)
+                    """{[el, idof, jdof]:
+                        0<=el<nelements and
+                        0<=idof, jdof<ndofs}
+                    """,
+                    """
+                    result[el, idof] = %s(jdof, operand[el, jdof])
+                    """ % reduction_name,
+                    default_offset=lp.auto,
+                    lang_version=MOST_RECENT_LANGUAGE_VERSION)
+
             knl = lp.tag_inames(knl, "el:g.0,idof:l.0")
             return knl
 
-        discr = self.bound_expr.get_discretization(expr.where)
+        @memoize_in(self.places, "elementwise_"+reduction_name)
+        def element_knl():
+            import loopy as lp
+            knl = lp.make_kernel(
+                    """{[el, jdof]:
+                        0<=el<nelements and
+                        0<=jdof<ndofs}
+                    """,
+                    """
+                    result[el] = %s(jdof, operand[el, jdof])
+                    """ % reduction_name,
+                    default_offset=lp.auto,
+                    lang_version=MOST_RECENT_LANGUAGE_VERSION)
 
-        operand = self.rec(expr.operand)
+            return knl
 
-        assert operand.shape == (discr.nnodes,)
+        def _reduce(nresult, knl, view):
+            result = cl.array.empty(self.queue, nresult, operand.dtype)
+            for igrp, group in enumerate(discr.groups):
+                knl(self.queue,
+                        operand=group.view(operand),
+                        result=view(igrp, result))
 
-        result = cl.array.empty(self.queue, discr.nnodes, operand.dtype)
-        for group in discr.groups:
-            knl()(self.queue,
-                    input=group.view(operand),
-                    result=group.view(result))
+            return result
 
-        return result
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
+        operand = self.rec(expr.operand)
+        assert operand.shape == (discr.nnodes,)
+
+        granularity = expr.dofdesc.granularity
+        if granularity is sym.GRANULARITY_NODE:
+            return _reduce(discr.nnodes,
+                    node_knl(),
+                    lambda g, x: discr.groups[g].view(x))
+        elif granularity is sym.GRANULARITY_ELEMENT:
+            return _reduce(discr.mesh.nelements,
+                    element_knl(),
+                    lambda g, x: mesh_el_view(discr.mesh, g, x))
+        else:
+            raise ValueError('unsupported granularity: %s' % granularity)
 
     def map_elementwise_sum(self, expr):
         return self._map_elementwise_reduction("sum", expr)
@@ -119,7 +169,7 @@ class EvaluationMapper(EvaluationMapperBase):
         return self._map_elementwise_reduction("max", expr)
 
     def map_ones(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
 
         result = (discr
                 .empty(queue=self.queue, dtype=discr.real_dtype)
@@ -129,12 +179,12 @@ class EvaluationMapper(EvaluationMapperBase):
         return result
 
     def map_node_coordinate_component(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
         return discr.nodes()[expr.ambient_axis] \
                 .with_queue(self.queue)
 
     def map_num_reference_derivative(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
 
         from pytools import flatten
         ref_axes = flatten([axis] * mult for axis, mult in expr.ref_axes)
@@ -144,7 +194,7 @@ class EvaluationMapper(EvaluationMapperBase):
                         .with_queue(self.queue)
 
     def map_q_weight(self, expr):
-        discr = self.bound_expr.get_discretization(expr.where)
+        discr = self.bound_expr.get_discretization(expr.dofdesc)
         return discr.quad_weights(self.queue) \
                 .with_queue(self.queue)
 
@@ -156,11 +206,11 @@ class EvaluationMapper(EvaluationMapperBase):
         except KeyError:
             bound_op = bind(
                     expr.expression,
-                    self.bound_expr.places[expr.where],
+                    self.places.get_geometry(expr.dofdesc),
                     self.bound_expr.iprec)
             bound_op_cache[expr] = bound_op
 
-        scipy_op = bound_op.scipy_op(expr.variable_name, expr.where,
+        scipy_op = bound_op.scipy_op(expr.variable_name, expr.dofdesc,
                 **dict((var_name, self.rec(var_expr))
                     for var_name, var_expr in six.iteritems(expr.extra_vars)))
 
@@ -170,35 +220,14 @@ class EvaluationMapper(EvaluationMapperBase):
         return result
 
     def map_interpolation(self, expr):
-        if not isinstance(expr.target, sym.QBXSourceQuadStage2):
-            raise RuntimeError("can only interpolate to quad_stage2 mesh")
-
-        source = self.bound_expr.places[expr.source]
         operand = self.rec(expr.operand)
 
-        from pytential.source import LayerPotentialSourceBase
-        if isinstance(source, LayerPotentialSourceBase):
-            sym_source = expr.source
-            if not isinstance(sym_source, sym._QBXSource):
-                sym_source = sym.QBXSourceStage1(sym_source)
-
-            if isinstance(sym_source, sym.QBXSourceStage1):
-                resampler = source.resampler
-            elif isinstance(sym_source, sym.QBXSourceStage2):
-                resampler = source.refined_interp_to_ovsmp_quad_connection
-            elif isinstance(sym_source, sym.QBXSourceQuadStage2):
-                resampler = lambda x, y: y  # noqa
-            else:
-                from pytential.symbolic.mappers import stringify_where
-                raise ValueError(
-                    "unknown `where` identifier in "
-                    "interpolation source: {}".format(
-                        stringify_where(sym_source)))
-        else:
-            raise TypeError("source must be a `LayerPotentialSourceBase`")
-
         if isinstance(operand, cl.array.Array):
-            return resampler(self.queue, operand)
+            from pytential.symbolic.dof_connection import connection_from_dds
+
+            conn = connection_from_dds(self.places,
+                    expr.from_dd, expr.to_dd)
+            return conn(self.queue, operand).with_queue(self.queue)
         elif isinstance(operand, (int, float, complex, np.number)):
             return operand
         else:
@@ -329,16 +358,15 @@ def _prepare_domains(nresults, places, domains, default_domain):
     """
 
     if domains is None:
-        if default_domain not in places:
-            raise RuntimeError("'domains is None' requires "
-                               "default domain to be defined in places")
         dom_name = default_domain
         return nresults * [dom_name]
     elif not isinstance(domains, (list, tuple)):
         dom_name = domains
         return nresults * [dom_name]
 
+    domains = [sym.as_dofdesc(d) for d in domains]
     assert len(domains) == nresults
+
     return domains
 
 
@@ -352,21 +380,26 @@ def _prepare_expr(places, expr):
 
     from pytential.source import LayerPotentialSourceBase
     from pytential.symbolic.mappers import (
-            ToTargetTagger, DerivativeBinder)
+            ToTargetTagger,
+            DerivativeBinder,
+            InterpolationPreprocessor)
 
-    expr = ToTargetTagger(*places._default_place_ids)(expr)
+    expr = ToTargetTagger(*places.auto_where)(expr)
     expr = DerivativeBinder()(expr)
 
     for name, place in six.iteritems(places.places):
         if isinstance(place, LayerPotentialSourceBase):
             expr = place.preprocess_optemplate(name, places, expr)
 
+    # NOTE: only insert interpolation operators after the layer potential
+    # operators were preprocessed to avoid any confusion
+    expr = InterpolationPreprocessor(places)(expr)
     return expr
 
 # }}}
 
 
-# {{{ bound expression
+# {{{ geometry collection
 
 class GeometryCollection(object):
     """A mapping from symbolic identifiers ("place IDs", typically strings)
@@ -378,8 +411,10 @@ class GeometryCollection(object):
     of subsets of them, as well as related common subexpressions such as
     metric terms.
 
-    .. method:: __getitem__
-    .. method:: get_discretization
+    .. automethod:: get_discretization
+    .. automethod:: get_geometry
+    .. automethod:: copy
+
     .. method:: get_cache
     """
 
@@ -400,31 +435,44 @@ class GeometryCollection(object):
         """
 
         from pytential.target import TargetBase
+        from pytential.source import PotentialSource
+        from pytential.qbx import QBXLayerPotentialSource
         from meshmode.discretization import Discretization
-        from pytential.source import LayerPotentialSourceBase, PotentialSource
 
-        if auto_where is None:
-            source_where, target_where = sym.DEFAULT_SOURCE, sym.DEFAULT_TARGET
+        # {{{ define default source and target descriptors
+
+        if isinstance(auto_where, (list, tuple)):
+            auto_source, auto_target = auto_where
         else:
-            # NOTE: keeping this here to make sure auto_where unpacks into
-            # just the two elements
-            source_where, target_where = auto_where
+            auto_source, auto_target = auto_where, None
+
+        if auto_source is None:
+            auto_source = sym.DEFAULT_SOURCE
+        if auto_target is None:
+            auto_target = sym.DEFAULT_TARGET
+
+        auto_source = sym.as_dofdesc(auto_source)
+        auto_target = sym.as_dofdesc(auto_target)
+        self.auto_where = (auto_source, auto_target)
+
+        # }}}
 
-        self._default_source_place = source_where
-        self._default_target_place = target_where
-        self._default_place_ids = (source_where, target_where)
+        # {{{ construct dict
 
         self.places = {}
-        if isinstance(places, LayerPotentialSourceBase):
-            self.places[source_where] = places
-            self.places[target_where] = \
-                    self._get_lpot_discretization(places, target_where)
-        elif isinstance(places, (Discretization, TargetBase)):
-            self.places[target_where] = places
+        if isinstance(places, QBXLayerPotentialSource):
+            self.places[auto_source.geometry] = places
+            self.places[auto_target.geometry] = \
+                    self._get_lpot_discretization(places, auto_target)
+        elif isinstance(places, (Discretization, PotentialSource)):
+            self.places[auto_source.geometry] = places
+            self.places[auto_target.geometry] = places
+        elif isinstance(places, TargetBase):
+            self.places[auto_target.geometry] = places
         elif isinstance(places, tuple):
             source_discr, target_discr = places
-            self.places[source_where] = source_discr
-            self.places[target_where] = target_discr
+            self.places[auto_source.geometry] = source_discr
+            self.places[auto_target.geometry] = target_discr
         else:
             self.places = places.copy()
 
@@ -433,62 +481,76 @@ class GeometryCollection(object):
                 raise TypeError("Must pass discretization, targets or "
                         "layer potential sources as 'places'.")
 
+        # }}}
+
         self.caches = {}
 
-    def _get_lpot_discretization(self, lpot, where):
-        if not isinstance(where, sym._QBXSource):
-            where = sym.QBXSourceStage1(where)
+    @property
+    def auto_source(self):
+        return self.auto_where[0]
 
-        if isinstance(where, sym.QBXSourceStage1):
-            return lpot.density_discr
-        if isinstance(where, sym.QBXSourceStage2):
+    @property
+    def auto_target(self):
+        return self.auto_where[1]
+
+    def _get_lpot_discretization(self, lpot, dofdesc):
+        if dofdesc.discr_stage == sym.QBX_SOURCE_STAGE2:
             return lpot.stage2_density_discr
-        if isinstance(where, sym.QBXSourceQuadStage2):
+        if dofdesc.discr_stage == sym.QBX_SOURCE_QUAD_STAGE2:
             return lpot.quad_stage2_density_discr
+        return lpot.density_discr
 
-        from pytential.symbolic.mappers import stringify_where
-        raise ValueError('unknown `where` identifier: {}'.format(
-            stringify_where(where)))
-
-    def get_discretization(self, where):
+    def get_discretization(self, dofdesc):
         """
-        :arg where: location identifier.
+        :arg dofdesc: a :class:`~pytential.symbolic.primitives.DOFDescriptor`
+            specifying the desired discretization.
 
         :return: a geometry object in the collection corresponding to the
-            key *where*. If it is a
+            key *dofdesc*. If it is a
             :class:`~pytential.source.LayerPotentialSourceBase`, we look for
             the corresponding :class:`~meshmode.discretization.Discretization`
             in its attributes instead.
         """
 
-        if where in self.places:
-            discr = self.places[where]
+        dofdesc = sym.as_dofdesc(dofdesc)
+        if dofdesc.geometry in self.places:
+            discr = self.places[dofdesc.geometry]
         else:
-            discr = self.places.get(getattr(where, 'where', None), None)
-
-        if discr is None:
-            raise KeyError('`where` not in the collection: {}'.format(where))
+            raise KeyError('geometry not in the collection: {}'.format(
+                dofdesc.geometry))
 
+        from pytential.qbx import QBXLayerPotentialSource
         from pytential.source import LayerPotentialSourceBase
-        if isinstance(discr, LayerPotentialSourceBase):
-            return self._get_lpot_discretization(discr, where)
+
+        if isinstance(discr, QBXLayerPotentialSource):
+            return self._get_lpot_discretization(discr, dofdesc)
+        elif isinstance(discr, LayerPotentialSourceBase):
+            return discr.density_discr
         else:
             return discr
 
-    def __getitem__(self, where):
-        return self.places[where]
-
-    def __contains__(self, where):
-        return where in self.places
+    def get_geometry(self, dofdesc):
+        dofdesc = sym.as_dofdesc(dofdesc)
+        return self.places[dofdesc.geometry]
 
     def copy(self):
         return GeometryCollection(
                 self.places,
-                auto_where=self._default_place_ids)
+                auto_where=self.auto_where)
 
     def get_cache(self, name):
         return self.caches.setdefault(name, {})
 
+    def __repr__(self):
+        return "%s(%s)" % (type(self).__name__, repr(self.places))
+
+    def __str__(self):
+        return "%s(%s)" % (type(self).__name__, str(self.places))
+
+# }}}
+
+
+# {{{ bound expression
 
 class BoundExpression(object):
     def __init__(self, places, sym_op_expr):
@@ -511,7 +573,7 @@ class BoundExpression(object):
             *None* values indicating the domains on which each component of the
             solution vector lives.  *None* values indicate that the component
             is a scalar.  If *domains* is *None*,
-            :class:`pytential.symbolic.primitives.DEFAULT_TARGET`, is required
+            :class:`~pytential.symbolic.primitives.DEFAULT_TARGET` is required
             to be a key in :attr:`places`.
         """
 
@@ -521,8 +583,8 @@ class BoundExpression(object):
         else:
             nresults = 1
 
-        domains = _prepare_domains(nresults, self.places, domains,
-                sym.DEFAULT_TARGET)
+        domains = _prepare_domains(nresults,
+                self.places, domains, self.places.auto_target)
 
         total_dofs = 0
         starts_and_ends = []
@@ -530,7 +592,7 @@ class BoundExpression(object):
             if dom_name is None:
                 size = 1
             else:
-                size = self.places[dom_name].nnodes
+                size = self.places.get_geometry(dom_name).nnodes
 
             starts_and_ends.append((total_dofs, total_dofs+size))
             total_dofs += size
@@ -559,10 +621,7 @@ def bind(places, expr, auto_where=None):
 
     if not isinstance(places, GeometryCollection):
         places = GeometryCollection(places, auto_where=auto_where)
-
-    from pytential.symbolic.mappers import QBXInterpolationPreprocessor
     expr = _prepare_expr(places, expr)
-    expr = QBXInterpolationPreprocessor(places)(expr)
 
     return BoundExpression(places, expr)
 
@@ -638,8 +697,8 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None,
         # not iterable, wrap in a list
         input_exprs = [input_exprs]
 
-    domains = _prepare_domains(len(input_exprs), places, domains,
-                               places._default_source_place)
+    domains = _prepare_domains(len(input_exprs),
+            places, domains, places.auto_source)
 
     from pytential.symbolic.matrix import MatrixBuilder, is_zero
     nblock_rows = len(exprs)
@@ -653,7 +712,7 @@ def build_matrix(queue, places, exprs, input_exprs, domains=None,
                 dep_expr=input_exprs[ibcol],
                 other_dep_exprs=(input_exprs[:ibcol]
                                  + input_exprs[ibcol + 1:]),
-                dep_source=places[domains[ibcol]],
+                dep_source=places.get_geometry(domains[ibcol]),
                 dep_discr=places.get_discretization(domains[ibcol]),
                 places=places,
                 context=context)
diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py
index 415453537e877bb40d7dbe35ade475e0337e51ef..66c0bba0f9f851c805357cf1f1d61af8c9b6415d 100644
--- a/pytential/symbolic/mappers.py
+++ b/pytential/symbolic/mappers.py
@@ -51,7 +51,6 @@ from pymbolic.geometric_algebra.mapper import (
         DerivativeSourceFinder
         as DerivativeSourceFinderBase,
 
-
         GraphvizMapper as GraphvizMapperBase)
 import pytential.symbolic.primitives as prim
 
@@ -63,14 +62,14 @@ class IdentityMapper(IdentityMapperBase):
     map_node_max = map_node_sum
 
     def map_elementwise_sum(self, expr):
-        return type(expr)(self.rec(expr.operand), expr.where)
+        return type(expr)(self.rec(expr.operand), expr.dofdesc)
 
     map_elementwise_min = map_elementwise_sum
     map_elementwise_max = map_elementwise_sum
 
     def map_num_reference_derivative(self, expr):
         return type(expr)(expr.ref_axes, self.rec(expr.operand),
-                expr.where)
+                expr.dofdesc)
 
     # {{{ childless -- no need to rebuild
 
@@ -93,7 +92,7 @@ class IdentityMapper(IdentityMapperBase):
                 dict([
                     (name, self.rec(name_expr))
                     for name, name_expr in six.iteritems(expr.extra_vars)]),
-                expr.where)
+                expr.dofdesc)
 
     def map_int_g(self, expr):
         return expr.copy(
@@ -104,7 +103,7 @@ class IdentityMapper(IdentityMapperBase):
                     ))
 
     def map_interpolation(self, expr):
-        return type(expr)(expr.source, expr.target, self.rec(expr.operand))
+        return type(expr)(expr.from_dd, expr.to_dd, self.rec(expr.operand))
 
 
 class CombineMapper(CombineMapperBase):
@@ -153,7 +152,7 @@ class DependencyMapper(DependencyMapperBase, Collector):
 class EvaluationMapper(EvaluationMapperBase):
     """Unlike :mod:`pymbolic.mapper.evaluation.EvaluationMapper`, this class
     does evaluation mostly to get :class:`pymbolic.geometric_algebra.MultiVector`
-    instances to to do their thing, and perhaps to automatically kill terms
+    instances to do their thing, and perhaps to automatically kill terms
     that are multiplied by zero. Otherwise it intends to largely preserve
     the structure of the input expression.
     """
@@ -178,7 +177,7 @@ class EvaluationMapper(EvaluationMapperBase):
     def map_num_reference_derivative(self, expr):
         return componentwise(
                 lambda subexpr: type(expr)(
-                        expr.ref_axes, self.rec(subexpr), expr.where),
+                        expr.ref_axes, self.rec(subexpr), expr.dofdesc),
                 expr.operand)
 
     def map_int_g(self, expr):
@@ -200,12 +199,11 @@ class EvaluationMapper(EvaluationMapperBase):
                 expr.scope)
 
 
-# {{{ target/source tagging
+# {{{ dofdesc tagging
 
 class LocationTagger(CSECachingMapperMixin, IdentityMapper):
-    """Used internally by :class:`ToTargetTagger`. Tags all src/target taggable
-    leaves/operators as going back onto the source.
-    """
+    """Used internally by :class:`ToTargetTagger`."""
+
     def __init__(self, default_where, default_source=prim.DEFAULT_SOURCE):
         self.default_source = default_source
         self.default_where = default_where
@@ -213,51 +211,54 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper):
     map_common_subexpression_uncached = \
             IdentityMapper.map_common_subexpression
 
+    def _default_dofdesc(self, dofdesc):
+        if dofdesc.geometry is None:
+            if dofdesc.discr_stage is None \
+                    and dofdesc.granularity == prim.GRANULARITY_NODE:
+                dofdesc = dofdesc.copy(geometry=self.default_where)
+            else:
+                dofdesc = dofdesc.copy(geometry=self.default_source)
+
+        return dofdesc
+
     def map_ones(self, expr):
-        if expr.where is None:
-            return type(expr)(where=self.default_where)
-        else:
-            return expr
+        return type(expr)(dofdesc=self._default_dofdesc(expr.dofdesc))
 
     map_q_weight = map_ones
 
     def map_parametrization_derivative_component(self, expr):
-        if expr.where is None:
-            return type(expr)(
-                    expr.ambient_axis, expr.ref_axis, self.default_where)
-        else:
-            return expr
+        return type(expr)(
+                expr.ambient_axis,
+                expr.ref_axis,
+                self._default_dofdesc(expr.dofdesc))
 
     def map_node_coordinate_component(self, expr):
-        if expr.where is None:
-            return type(expr)(
-                    expr.ambient_axis, self.default_where)
-        else:
-            return expr
+        return type(expr)(
+                expr.ambient_axis,
+                self._default_dofdesc(expr.dofdesc))
 
     def map_num_reference_derivative(self, expr):
-        if expr.where is None:
-            return type(expr)(
-                    expr.ref_axes, self.rec(expr.operand), self.default_where)
-        else:
-            return expr
+        return type(expr)(
+                expr.ref_axes,
+                self.rec(expr.operand),
+                self._default_dofdesc(expr.dofdesc))
 
     def map_elementwise_sum(self, expr):
         return type(expr)(
                 self.rec(expr.operand),
-                expr.where if expr.where is not None else self.default_where)
+                self._default_dofdesc(expr.dofdesc))
 
     map_elementwise_min = map_elementwise_sum
     map_elementwise_max = map_elementwise_sum
 
     def map_int_g(self, expr):
         source = expr.source
-        target = expr.target
+        if source.geometry is None:
+            source = source.copy(geometry=self.default_source)
 
-        if source is None:
-            source = self.default_source
-        if target is None:
-            target = self.default_where
+        target = expr.target
+        if target.geometry is None:
+            target = target.copy(geometry=self.default_where)
 
         return type(expr)(
                 expr.kernel,
@@ -269,10 +270,9 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper):
                     ))
 
     def map_inverse(self, expr):
-        where = expr.where
-
-        if where is None:
-            where = self.default_where
+        dofdesc = expr.dofdesc
+        if dofdesc.geometry is None:
+            dofdesc = dofdesc.copy(geometry=self.default_where)
 
         return type(expr)(
                 # don't recurse into expression--it is a separate world that
@@ -282,28 +282,35 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper):
                 dict([
                     (name, self.rec(name_expr))
                     for name, name_expr in six.iteritems(expr.extra_vars)]),
-                where)
+                dofdesc)
 
     def map_interpolation(self, expr):
-        source = expr.source
-        target = expr.target
+        from_dd = expr.from_dd
+        if from_dd.geometry is None:
+            from_dd = from_dd.copy(geometry=self.default_source)
 
-        if source is None:
-            source = self.default_source
-        if target is None:
-            target = self.default_where
+        to_dd = expr.to_dd
+        if to_dd.geometry is None:
+            to_dd = to_dd.copy(geometry=self.default_source)
 
-        return type(expr)(source, target, expr.operand)
+        return type(expr)(from_dd, to_dd, self.operand_rec(expr.operand))
 
     def operand_rec(self, expr):
         return self.rec(expr)
 
 
 class ToTargetTagger(LocationTagger):
-    """Descends into the expression tree, marking everything up to the first
-    layer potential operator as operating on the targets, and everything below
-    there as operating on the source. Also performs a consistency check such
-    that only 'target' and 'source' items are combined arithmetically.
+    """Descends into the expression tree, marking expressions based on two
+    heuristics:
+
+    * everything up to the first layer potential operator is marked as
+    operating on the targets, and everything below there as operating on the
+    source.
+    * if an expression has a :class:`~pytential.symbolic.primitives.DOFDescriptor`
+    that requires a :class:`~pytential.source.LayerPotentialSourceBase` to be
+    used (e.g. by being defined on
+    :class:`~pytential.symbolic.primitives.QBX_SOURCE_QUAD_STAGE2`), then
+    it is marked as operating on a source.
     """
 
     def __init__(self, default_source, default_target):
@@ -312,6 +319,46 @@ class ToTargetTagger(LocationTagger):
         self.operand_rec = LocationTagger(default_source,
                                           default_source=default_source)
 
+
+class DiscretizationStageTagger(IdentityMapper):
+    """Descends into an expression tree and changes the
+    :attr:`~pytential.symbolic.primitives.DOFDescriptor.discr_stage` to
+    :attr:`discr_stage`.
+
+    .. attribute:: discr_stage
+
+        The new discretization for the DOFs in the expression. For valid
+        values, see
+        :attr:`~pytential.symbolic.primitives.DOFDescriptor.discr_stage`.
+    """
+
+    def __init__(self, discr_stage):
+        if not (discr_stage == prim.QBX_SOURCE_STAGE1
+                or discr_stage == prim.QBX_SOURCE_STAGE2
+                or discr_stage == prim.QBX_SOURCE_QUAD_STAGE2):
+            raise ValueError('unknown discr stage tag: "{}"'.format(discr_stage))
+
+        self.discr_stage = discr_stage
+
+    def map_node_coordinate_component(self, expr):
+        dofdesc = expr.dofdesc
+        if dofdesc.discr_stage == self.discr_stage:
+            return expr
+
+        return type(expr)(
+                expr.ambient_axis,
+                dofdesc.copy(discr_stage=self.discr_stage))
+
+    def map_num_reference_derivative(self, expr):
+        dofdesc = expr.dofdesc
+        if dofdesc.discr_stage == self.discr_stage:
+            return expr
+
+        return type(expr)(
+                expr.ref_axes,
+                self.rec(expr.operand),
+                dofdesc.copy(discr_stage=self.discr_stage))
+
 # }}}
 
 
@@ -406,36 +453,66 @@ class UnregularizedPreprocessor(IdentityMapper):
 # }}}
 
 
-# {{{ QBX preprocessor
+# {{{ interpolation preprocessor
+
+class InterpolationPreprocessor(IdentityMapper):
+    """Handle expressions that require upsampling or downsampling by inserting
+    a :class:`~pytential.symbolic.primitives.Interpolation`. This is used to
+
+    * do differentiation on
+    :attr:`~pytential.source.LayerPotentialSource.quad_stage2_density_discr`,
+    by performing it on
+    :attr:`~pytential.source.LayerPotentialSource.stage2_density_discr` and
+    upsampling.
+    * upsample layer potential sources to
+    :attr:`~pytential.source.LayerPotentialSource.quad_stage2_density_discr`,
+    """
 
-class QBXInterpolationPreprocessor(IdentityMapper):
     def __init__(self, places):
         self.places = places
+        self.from_discr_stage = prim.QBX_SOURCE_STAGE2
+        self.tagger = DiscretizationStageTagger(self.from_discr_stage)
+
+    def map_num_reference_derivative(self, expr):
+        to_dd = expr.dofdesc
+        if to_dd.discr_stage != prim.QBX_SOURCE_QUAD_STAGE2:
+            return expr
+
+        from pytential.qbx import QBXLayerPotentialSource
+        lpot_source = self.places.get_geometry(to_dd)
+        if not isinstance(lpot_source, QBXLayerPotentialSource):
+            return expr
+
+        from_dd = to_dd.copy(discr_stage=self.from_discr_stage)
+        return prim.interp(from_dd, to_dd, self.rec(self.tagger(expr)))
 
     def map_int_g(self, expr):
-        if isinstance(expr.source, prim.QBXSourceQuadStage2):
+        from_dd = expr.source
+        if from_dd.discr_stage is not None:
             return expr
 
-        from pytential.source import LayerPotentialSourceBase
-        source = self.places[expr.source]
+        from pytential.qbx import QBXLayerPotentialSource
+        lpot_source = self.places.get_geometry(from_dd)
+        if not isinstance(lpot_source, QBXLayerPotentialSource):
+            return expr
 
-        if isinstance(source, LayerPotentialSourceBase):
-            stage2_source = prim.QBXSourceQuadStage2(prim.DEFAULT_SOURCE)
+        to_dd = from_dd.copy(discr_stage=prim.QBX_SOURCE_QUAD_STAGE2)
+        density = prim.interp(from_dd, to_dd, self.rec(expr.density))
+        kernel_arguments = dict(
+                (name, prim.interp(from_dd, to_dd, self.rec(arg_expr)))
+                for name, arg_expr in expr.kernel_arguments.items())
 
-            density = prim.Interpolation(
-                    expr.source, stage2_source, self.rec(expr.density))
-            kernel_arguments = dict(
-                    (name, prim.Interpolation(
-                        expr.source, stage2_source, self.rec(arg_expr)))
-                    for name, arg_expr in expr.kernel_arguments.items())
+        return expr.copy(
+                kernel=expr.kernel,
+                density=density,
+                kernel_arguments=kernel_arguments,
+                source=to_dd,
+                target=expr.target)
 
-            expr = expr.copy(
-                    kernel=expr.kernel,
-                    density=density,
-                    kernel_arguments=kernel_arguments)
+# }}}
 
-        return expr
 
+# {{{ QBX preprocessor
 
 class QBXPreprocessor(IdentityMapper):
     def __init__(self, source_name, places):
@@ -443,7 +520,10 @@ class QBXPreprocessor(IdentityMapper):
         self.places = places
 
     def map_int_g(self, expr):
-        source_discr = self.places.get_discretization(self.source_name)
+        if expr.source.geometry != self.source_name:
+            return expr
+
+        source_discr = self.places.get_discretization(expr.source)
         target_discr = self.places.get_discretization(expr.target)
 
         if expr.qbx_forced_limit == 0:
@@ -495,27 +575,13 @@ class QBXPreprocessor(IdentityMapper):
 # {{{ stringifier
 
 def stringify_where(where):
-    if isinstance(where, prim.QBXSourceStage1):
-        return "stage1(%s)" % stringify_where(where.where)
-    if isinstance(where, prim.QBXSourceStage2):
-        return "stage2(%s)" % stringify_where(where.where)
-    if isinstance(where, prim.QBXSourceQuadStage2):
-        return "quad_stage2(%s)" % stringify_where(where.where)
-
-    if where is None:
-        return "?"
-    elif where is prim.DEFAULT_SOURCE:
-        return "s"
-    elif where is prim.DEFAULT_TARGET:
-        return "t"
-    else:
-        return str(where)
+    return str(prim.as_dofdesc(where))
 
 
 class StringifyMapper(BaseStringifyMapper):
 
     def map_ones(self, expr, enclosing_prec):
-        return "Ones.%s" % stringify_where(expr.where)
+        return "Ones[%s]" % stringify_where(expr.dofdesc)
 
     def map_inverse(self, expr, enclosing_prec):
         return "Solve(%s = %s {%s})" % (
@@ -532,18 +598,18 @@ class StringifyMapper(BaseStringifyMapper):
                 set())
 
     def map_elementwise_sum(self, expr, enclosing_prec):
-        return "ElwiseSum.%s(%s)" % (
-                stringify_where(expr.where),
+        return "ElwiseSum[%s](%s)" % (
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_NONE))
 
     def map_elementwise_min(self, expr, enclosing_prec):
-        return "ElwiseMin.%s(%s)" % (
-                stringify_where(expr.where),
+        return "ElwiseMin[%s](%s)" % (
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_NONE))
 
     def map_elementwise_max(self, expr, enclosing_prec):
-        return "ElwiseMax.%s(%s)" % (
-                stringify_where(expr.where),
+        return "ElwiseMax[%s](%s)" % (
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_NONE))
 
     def map_node_max(self, expr, enclosing_prec):
@@ -553,8 +619,8 @@ class StringifyMapper(BaseStringifyMapper):
         return "NodeSum(%s)" % self.rec(expr.operand, PREC_NONE)
 
     def map_node_coordinate_component(self, expr, enclosing_prec):
-        return "x%d.%s" % (expr.ambient_axis,
-                stringify_where(expr.where))
+        return "x%d[%s]" % (expr.ambient_axis,
+                stringify_where(expr.dofdesc))
 
     def map_num_reference_derivative(self, expr, enclosing_prec):
         diff_op = " ".join(
@@ -563,9 +629,9 @@ class StringifyMapper(BaseStringifyMapper):
                 "d/dr%d^%d" % (axis, mult)
                 for axis, mult in expr.ref_axes)
 
-        result = "%s.%s %s" % (
+        result = "%s[%s] %s" % (
                 diff_op,
-                stringify_where(expr.where),
+                stringify_where(expr.dofdesc),
                 self.rec(expr.operand, PREC_PRODUCT),
                 )
 
@@ -575,10 +641,10 @@ class StringifyMapper(BaseStringifyMapper):
             return result
 
     def map_parametrization_derivative(self, expr, enclosing_prec):
-        return "dx/dr.%s" % (stringify_where(expr.where))
+        return "dx/dr[%s]" % (stringify_where(expr.dofdesc))
 
     def map_q_weight(self, expr, enclosing_prec):
-        return "w_quad.%s" % stringify_where(expr.where)
+        return "w_quad[%s]" % stringify_where(expr.dofdesc)
 
     def _stringify_kernel_args(self, kernel_arguments):
         if not kernel_arguments:
@@ -600,8 +666,8 @@ class StringifyMapper(BaseStringifyMapper):
 
     def map_interpolation(self, expr, enclosing_prec):
         return "Interp[%s->%s](%s)" % (
-                stringify_where(expr.source),
-                stringify_where(expr.target),
+                stringify_where(expr.from_dd),
+                stringify_where(expr.to_dd),
                 self.rec(expr.operand, PREC_PRODUCT))
 
 # }}}
diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py
index 9d0c2b2f66a7080f91cda5c668377c9985bfca3f..9f9e0e4a268b7f365003b9c293b7b08e246bf01c 100644
--- a/pytential/symbolic/matrix.py
+++ b/pytential/symbolic/matrix.py
@@ -43,34 +43,6 @@ def is_zero(x):
     return isinstance(x, (int, float, complex, np.number)) and x == 0
 
 
-def _resample_arg(queue, source, x):
-    """
-    :arg queue: a :class:`pyopencl.CommandQueue`.
-    :arg source: a :class:`pytential.source.LayerPotentialSourceBase` subclass.
-        If it is not a layer potential source, no resampling is done.
-    :arg x: a :class:`numpy.ndarray`.
-
-    :return: a resampled :class:`numpy.ndarray` (see
-        :method:`pytential.source.LayerPotentialSourceBase.resampler`).
-    """
-
-    from pytential.source import LayerPotentialSourceBase
-    if not isinstance(source, LayerPotentialSourceBase):
-        return x
-
-    if not isinstance(x, np.ndarray):
-        return x
-
-    if len(x.shape) >= 2:
-        raise RuntimeError("matrix variables in kernel arguments")
-
-    def resample(y):
-        return source.resampler(queue, cl.array.to_device(queue, y)).get(queue)
-
-    from pytools.obj_array import with_object_array_or_scalar
-    return with_object_array_or_scalar(resample, x)
-
-
 def _get_layer_potential_args(mapper, expr, source):
     """
     :arg mapper: a :class:`pytential.symbolic.matrix.MatrixBuilderBase`.
@@ -80,17 +52,9 @@ def _get_layer_potential_args(mapper, expr, source):
     :return: a mapping of kernel arguments evaluated by the *mapper*.
     """
 
-    # skip resampling if source and target are the same
-    from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET
-    if ((expr.source is not DEFAULT_SOURCE)
-            and (expr.target is not DEFAULT_TARGET)
-            and (isinstance(expr.source, type(expr.target)))):
-        source = None
-
     kernel_args = {}
     for arg_name, arg_expr in six.iteritems(expr.kernel_arguments):
-        rec_arg = mapper.rec(arg_expr)
-        kernel_args[arg_name] = _resample_arg(mapper.queue, source, rec_arg)
+        kernel_args[arg_name] = mapper.rec(arg_expr)
 
     return kernel_args
 
@@ -113,9 +77,7 @@ def _get_kernel_args(mapper, kernel, expr, source):
     for arg_name, arg_expr in six.iteritems(expr.kernel_arguments):
         if arg_name not in inner_kernel_args:
             continue
-
-        rec_arg = mapper.rec(arg_expr)
-        kernel_args[arg_name] = _resample_arg(mapper.queue, source, rec_arg)
+        kernel_args[arg_name] = mapper.rec(arg_expr)
 
     return kernel_args
 
@@ -154,13 +116,16 @@ def _get_centers_and_expansion_radii(queue, source, target_discr, qbx_forced_lim
 
     if source.density_discr is target_discr:
         # NOTE: skip expensive target association
-        from pytential.qbx.utils import get_centers_on_side
-        centers = get_centers_on_side(source, qbx_forced_limit)
-        radii = source._expansion_radii('nsources')
+        centers = bind(source,
+            sym.expansion_centers(source.ambient_dim, qbx_forced_limit))(queue)
+        radii = bind(source,
+            sym.expansion_radii(source.ambient_dim))(queue)
     else:
         from pytential.qbx.utils import get_interleaved_centers
         centers = get_interleaved_centers(queue, source)
-        radii = source._expansion_radii('ncenters')
+        radii = bind(source, sym.expansion_radii(
+            source.ambient_dim,
+            granularity=sym.GRANULARITY_CENTER))(queue)
 
         # NOTE: using a very small tolerance to make sure all the stage2
         # targets are associated to a center. We can't use the user provided
@@ -319,15 +284,16 @@ class MatrixBuilderBase(EvaluationMapperBase):
         if self.is_kind_matrix(rec_operand):
             raise NotImplementedError("derivatives")
 
-        where_discr = self.places[expr.where]
-        op = sym.NumReferenceDerivative(expr.ref_axes, sym.var("u"))
-        return bind(where_discr, op)(
-                self.queue, u=cl.array.to_device(self.queue, rec_operand)).get()
+        rec_operand = cl.array.to_device(self.queue, rec_operand)
+        op = sym.NumReferenceDerivative(
+                ref_axes=expr.ref_axes,
+                operand=sym.var("u"),
+                dofdesc=expr.dofdesc)
+        return bind(self.places, op)(self.queue, u=rec_operand).get()
 
     def map_node_coordinate_component(self, expr):
-        where_discr = self.places[expr.where]
-        op = sym.NodeCoordinateComponent(expr.ambient_axis)
-        return bind(where_discr, op)(self.queue).get()
+        op = sym.NodeCoordinateComponent(expr.ambient_axis, dofdesc=expr.dofdesc)
+        return bind(self.places, op)(self.queue).get()
 
     def map_call(self, expr):
         arg, = expr.parameters
@@ -401,13 +367,30 @@ class MatrixBuilder(MatrixBuilderBase):
         super(MatrixBuilder, self).__init__(queue, dep_expr, other_dep_exprs,
                 dep_source, dep_discr, places, context)
 
-    def map_int_g(self, expr):
-        where_source = expr.source
-        if where_source is sym.DEFAULT_SOURCE:
-            where_source = sym.QBXSourceQuadStage2(expr.source)
+    def map_interpolation(self, expr):
+        if expr.to_dd.discr_stage != sym.QBX_SOURCE_QUAD_STAGE2:
+            raise RuntimeError("can only interpolate to QBX_SOURCE_QUAD_STAGE2")
+
+        operand = self.rec(expr.operand)
+        if isinstance(operand, (int, float, complex, np.number)):
+            return operand
+        elif isinstance(operand, np.ndarray) and operand.ndim == 1:
+            from pytential.symbolic.dof_connection import connection_from_dds
+            conn = connection_from_dds(self.places,
+                    expr.from_dd, expr.to_dd)
+
+            operand = cl.array.to_device(self.queue, operand)
+            return conn(self.queue, operand).get(self.queue)
+        elif isinstance(operand, np.ndarray) and operand.ndim == 2:
+            resampler = self.places.get_geometry(expr.from_dd).direct_resampler
+            mat = resampler.full_resample_matrix(self.queue).get(self.queue)
+            return mat.dot(operand)
+        else:
+            raise RuntimeError('unknown operand type: {}'.format(type(operand)))
 
-        source = self.places[expr.source]
-        source_discr = self.places.get_discretization(where_source)
+    def map_int_g(self, expr):
+        lpot_source = self.places.get_geometry(expr.source)
+        source_discr = self.places.get_discretization(expr.source)
         target_discr = self.places.get_discretization(expr.target)
 
         rec_density = self.rec(expr.density)
@@ -419,10 +402,10 @@ class MatrixBuilder(MatrixBuilderBase):
             raise NotImplementedError("layer potentials on non-variables")
 
         kernel = expr.kernel
-        kernel_args = _get_layer_potential_args(self, expr, source)
+        kernel_args = _get_layer_potential_args(self, expr, lpot_source)
 
         from sumpy.expansion.local import LineTaylorLocalExpansion
-        local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order)
+        local_expn = LineTaylorLocalExpansion(kernel, lpot_source.qbx_order)
 
         from sumpy.qbx import LayerPotentialMatrixGenerator
         mat_gen = LayerPotentialMatrixGenerator(
@@ -430,7 +413,7 @@ class MatrixBuilder(MatrixBuilderBase):
 
         assert abs(expr.qbx_forced_limit) > 0
         centers, radii = _get_centers_and_expansion_radii(self.queue,
-                source, target_discr, expr.qbx_forced_limit)
+                lpot_source, target_discr, expr.qbx_forced_limit)
 
         _, (mat,) = mat_gen(self.queue,
                 targets=target_discr.nodes(),
@@ -440,17 +423,8 @@ class MatrixBuilder(MatrixBuilderBase):
                 **kernel_args)
         mat = mat.get()
 
-        waa = _get_weights_and_area_elements(self.queue, source, source_discr)
+        waa = _get_weights_and_area_elements(self.queue, lpot_source, source_discr)
         mat[:, :] *= waa.get(self.queue)
-
-        if target_discr.nnodes != source_discr.nnodes:
-            # NOTE: we only resample sources
-            assert target_discr.nnodes < source_discr.nnodes
-
-            resampler = source.direct_resampler
-            resample_mat = resampler.full_resample_matrix(self.queue).get(self.queue)
-            mat = mat.dot(resample_mat)
-
         mat = mat.dot(rec_density)
 
         return mat
@@ -470,7 +444,7 @@ class P2PMatrixBuilder(MatrixBuilderBase):
         self.exclude_self = exclude_self
 
     def map_int_g(self, expr):
-        source = self.places[expr.source]
+        source = self.places.get_geometry(expr.source)
         source_discr = self.places.get_discretization(expr.source)
         target_discr = self.places.get_discretization(expr.target)
 
@@ -533,7 +507,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase):
         return np.equal(tgtindices, srcindices).astype(np.float64)
 
     def map_int_g(self, expr):
-        source = self.places[expr.source]
+        source = self.places.get_geometry(expr.source)
         source_discr = self.places.get_discretization(expr.source)
         target_discr = self.places.get_discretization(expr.target)
 
@@ -548,7 +522,7 @@ class NearFieldBlockBuilder(MatrixBlockBuilderBase):
             raise NotImplementedError()
 
         kernel = expr.kernel
-        kernel_args = _get_layer_potential_args(self.mat_mapper, expr, source)
+        kernel_args = _get_layer_potential_args(self.mat_mapper, expr, None)
 
         from sumpy.expansion.local import LineTaylorLocalExpansion
         local_expn = LineTaylorLocalExpansion(kernel, source.qbx_order)
@@ -599,7 +573,7 @@ class FarFieldBlockBuilder(MatrixBlockBuilderBase):
         return np.equal(tgtindices, srcindices).astype(np.float64)
 
     def map_int_g(self, expr):
-        source = self.places[expr.source]
+        source = self.places.get_geometry(expr.source)
         source_discr = self.places.get_discretization(expr.source)
         target_discr = self.places.get_discretization(expr.target)
 
diff --git a/pytential/symbolic/pde/maxwell/__init__.py b/pytential/symbolic/pde/maxwell/__init__.py
index a8dc0c42390e7d56dc4220a7f13fdfe43e440264..d89d393c4948a2944d4efa4da2b89db65e55d4b1 100644
--- a/pytential/symbolic/pde/maxwell/__init__.py
+++ b/pytential/symbolic/pde/maxwell/__init__.py
@@ -70,7 +70,8 @@ def get_sym_maxwell_point_source(kernel, jxyz, k):
 
 # {{{ plane wave
 
-def get_sym_maxwell_plane_wave(amplitude_vec, v, omega, epsilon=1, mu=1, where=None):
+def get_sym_maxwell_plane_wave(amplitude_vec, v, omega,
+        epsilon=1, mu=1, dofdesc=None):
     r"""Return a symbolic expression that, when bound to a
     :class:`pytential.source.PointPotentialSource` will yield
     a field satisfying Maxwell's equations.
@@ -96,7 +97,7 @@ def get_sym_maxwell_plane_wave(amplitude_vec, v, omega, epsilon=1, mu=1, where=N
 
     # NOTE: for complex, need to ensure real(n).dot(imag(n)) = 0  (7.15)
 
-    x = sym.nodes(3, where).as_vector()
+    x = sym.nodes(3, dofdesc=dofdesc).as_vector()
 
     v_mag_squared = sym.cse(np.dot(v, v), "v_mag_squared")
     n = v/sym.sqrt(v_mag_squared)
diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py
index 840d804429344f16c62e8f388c654c65dffc6f5b..882b1b3c53cb53830bb94796b887a71494fdae32 100644
--- a/pytential/symbolic/pde/scalar.py
+++ b/pytential/symbolic/pde/scalar.py
@@ -48,15 +48,16 @@ class L2WeightedPDEOperator(object):
             warn("should use L2 weighting in %s" % type(self).__name__,
                     stacklevel=3)
 
-    def get_weight(self, where=None):
+    def get_weight(self, dofdesc=None):
         if self.use_l2_weighting:
-            return cse(area_element(self.kernel.dim, where=where)*QWeight(where))
+            return cse(area_element(self.kernel.dim, dofdesc=dofdesc)
+                    * QWeight(dofdesc=dofdesc))
         else:
             return 1
 
-    def get_sqrt_weight(self, where=None):
+    def get_sqrt_weight(self, dofdesc=None):
         if self.use_l2_weighting:
-            return sqrt_jac_q_weight(self.kernel.dim, where=where)
+            return sqrt_jac_q_weight(self.kernel.dim, dofdesc=dofdesc)
         else:
             return 1
 
diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py
index ec129369422f1a04d3f16144df47240b739269b7..1f4fd48583468dc31384429796527564cd424669 100644
--- a/pytential/symbolic/primitives.py
+++ b/pytential/symbolic/primitives.py
@@ -25,6 +25,7 @@ THE SOFTWARE.
 import six
 from six.moves import intern
 from warnings import warn
+from functools import wraps
 
 import numpy as np
 from pymbolic.primitives import (  # noqa: F401,N813
@@ -41,8 +42,11 @@ from functools import partial
 
 
 __doc__ = """
-.. |where-blurb| replace:: A symbolic name for a
-    :class:`pytential.discretization.Discretization`
+.. |dofdesc-blurb| replace:: A
+    :class:`~pytential.symbolic.primitives.DOFDescriptor` or a symbolic
+    name for a geometric object (such as a
+    :class:`~meshmode.discretization.Discretization`).
+
 
 Object types
 ^^^^^^^^^^^^
@@ -72,6 +76,23 @@ objects occur as part of a symbolic operator representation:
 those), which is not visible as an array axis in symbolic code. (They're
 visible only once evaluated.)
 
+DOF Description
+^^^^^^^^^^^^^^^
+
+.. autoclass:: DEFAULT_SOURCE
+.. autoclass:: DEFAULT_TARGET
+
+.. autoclass:: QBX_SOURCE_STAGE1
+.. autoclass:: QBX_SOURCE_STAGE2
+.. autoclass:: QBX_SOURCE_QUAD_STAGE2
+
+.. autoclass:: GRANULARITY_NODE
+.. autoclass:: GRANULARITY_CENTER
+.. autoclass:: GRANULARITY_ELEMENT
+
+.. autoclass:: DOFDescriptor
+.. autofunction:: as_dofdesc
+
 Placeholders
 ^^^^^^^^^^^^
 
@@ -126,6 +147,11 @@ Discretization properties
 .. autofunction:: second_fundamental_form
 .. autofunction:: shape_operator
 
+.. autofunction:: expansion_radii
+.. autofunction:: expansion_centers
+.. autofunction:: h_max
+.. autofunction:: weights_and_area_elements
+
 Elementary numerics
 ^^^^^^^^^^^^^^^^^^^
 
@@ -145,6 +171,7 @@ Operators
 ^^^^^^^^^
 
 .. autoclass:: Interpolation
+.. autofunction:: interp
 
 Geometric Calculus (based on Geometric/Clifford Algebra)
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -196,69 +223,198 @@ Pretty-printing expressions
 """
 
 
-# {{{ 'where' specifiers
+def _deprecate_kwargs(oldkey, newkey):
+    def super_wrapper(func):
+        @wraps(func)
+        def wrapper(*args, **kwargs):
+            if oldkey in kwargs:
+                warn("using `{}` keyword is deprecated. "
+                        "use `{}` instead".format(oldkey, newkey),
+                        DeprecationWarning, stacklevel=2)
+                kwargs[newkey] = kwargs[oldkey]
+                del kwargs[oldkey]
+
+            return func(*args, **kwargs)
+        return wrapper
 
-class DEFAULT_SOURCE:  # noqa
+    return super_wrapper
+
+
+# {{{ dof descriptors
+
+class DEFAULT_SOURCE:  # noqa: N801
     pass
 
 
-class DEFAULT_TARGET:  # noqa
+class DEFAULT_TARGET:  # noqa: N801
     pass
 
 
-class _QBXSource(object):
-    """A symbolic 'where' specifier for the a density of a
-    :attr:`pytential.qbx.QBXLayerPotentialSource`
-    layer potential source identified by :attr:`where`.
+class QBX_SOURCE_STAGE1:   # noqa: N801
+    """Symbolic identifier for the base `stage1` discretization
+    :attr:`pytential.source.LayerPotentialSourceBase.density_discr`.
+    """
+    pass
 
-    .. attribute:: where
 
-        An identifier of a layer potential source, as used in
-        :func:`pytential.bind`.
+class QBX_SOURCE_STAGE2:   # noqa: N801
+    """Symbolic identifier for the `stage2` discretization
+    :attr:`pytential.source.LayerPotentialSourceBase.stage2_density_discr`.
+    """
+    pass
 
-    .. note::
 
-        This is not documented functionality and only intended for
-        internal use.
+class QBX_SOURCE_QUAD_STAGE2:   # noqa: N801
+    """Symbolic identifier for the `stage2` discretization
+    :attr:`pytential.source.LayerPotentialSourceBase.quad_stage2_density_discr`.
     """
+    pass
+
+
+class GRANULARITY_NODE:     # noqa: N801
+    """DOFs are per node."""
+    pass
+
+
+class GRANULARITY_CENTER:   # noqa: N801
+    """DOFs interleaved per expansion center (two per node, one on each side)."""
+    pass
+
+
+class GRANULARITY_ELEMENT:  # noqa: N801
+    """DOFs per discretization element."""
+    pass
+
+
+class DOFDescriptor(object):
+    """A data structure specifying the meaning of a vector of degrees of freedom
+    that is handled by :mod:`pytential` (a "DOF vector"). In particular, using
+    :attr:`geometry`, this data structure describes the geometric object on which
+    the (scalar) function described by the DOF vector exists. Using
+    :attr:`granularity`, the data structure describes how the geometric object
+    is discretized (e.g. conventional nodal data, per-element scalars, etc.)
+
+    .. attribute:: geometry
+
+        An identifier for the geometry on which the DOFs exist. This can be a
+        simple string or any other hashable identifier for the geometric object.
+        The geometric objects are generally subclasses of
+        :class:`~pytential.source.PotentialSource`,
+        :class:`~pytential.target.TargetBase` or
+        :class:`~meshmode.discretization.Discretization`.
+
+    .. attribute:: discr_stage
+
+        Specific to a :class:`pytential.source.LayerPotentialSourceBase`,
+        this describes on which of the discretizations the
+        DOFs are defined. Can be one of :class:`QBX_SOURCE_STAGE1`,
+        :class:`QBX_SOURCE_STAGE2` or :class:`QBX_SOURCE_QUAD_STAGE2`.
 
-    def __init__(self, where):
-        self.where = where
+    .. attribute:: granularity
+
+        Describes the level of granularity of the DOF vector.
+        Can be one of :class:`GRANULARITY_NODE` (one DOF per node),
+        :class:`GRANULARITY_CENTER` (two DOFs per node, one per side) or
+        :class:`GRANULARITY_ELEMENT` (one DOF per element).
+    """
+
+    def __init__(self, geometry=None, discr_stage=None, granularity=None):
+        if granularity is None:
+            granularity = GRANULARITY_NODE
+
+        if discr_stage is not None:
+            if not (discr_stage == QBX_SOURCE_STAGE1
+                    or discr_stage == QBX_SOURCE_STAGE2
+                    or discr_stage == QBX_SOURCE_QUAD_STAGE2):
+                raise ValueError('unknown discr stage tag: "{}"'.format(discr_stage))
+
+        if not (granularity == GRANULARITY_NODE
+                or granularity == GRANULARITY_CENTER
+                or granularity == GRANULARITY_ELEMENT):
+            raise ValueError('unknown granularity: "{}"'.format(granularity))
+
+        self.geometry = geometry
+        self.discr_stage = discr_stage
+        self.granularity = granularity
+
+    def copy(self, geometry=None, discr_stage=None, granularity=None):
+        if isinstance(geometry, DOFDescriptor):
+            discr_stage = geometry.discr_stage \
+                    if discr_stage is None else discr_stage
+            geometry = geometry.geometry
+
+        return type(self)(
+                geometry=(self.geometry
+                    if geometry is None else geometry),
+                granularity=(self.granularity
+                    if granularity is None else granularity),
+                discr_stage=(self.discr_stage
+                    if discr_stage is None else discr_stage))
 
     def __hash__(self):
-        return hash((type(self), self.where))
+        return hash((type(self),
+            self.geometry, self.discr_stage, self.granularity))
 
     def __eq__(self, other):
-        return type(self) is type(other) and self.where == other.where
+        return (type(self) is type(other)
+                and self.geometry == other.geometry
+                and self.discr_stage == other.discr_stage
+                and self.granularity == other.granularity)
 
     def __ne__(self, other):
         return not self.__eq__(other)
 
+    def __repr__(self):
+        discr_stage = self.discr_stage \
+                if self.discr_stage is None else self.discr_stage.__name__,
+        granularity = self.granularity.__name__
+        return '{}(geometry={}, stage={}, granularity={})'.format(
+                type(self).__name__, self.geometry, discr_stage, granularity)
+
+    def __str__(self):
+        name = []
+        if self.geometry is None:
+            name.append("?")
+        elif self.geometry == DEFAULT_SOURCE:
+            name.append("s")
+        elif self.geometry == DEFAULT_TARGET:
+            name.append("t")
+        else:
+            name.append(str(self.geometry))
+
+        if self.discr_stage == QBX_SOURCE_STAGE2:
+            name.append("stage2")
+        elif self.discr_stage == QBX_SOURCE_QUAD_STAGE2:
+            name.append("quads2")
 
-class QBXSourceStage1(_QBXSource):
-    """An explicit symbolic 'where' specifier for the
-    :attr:`pytential.qbx.QBXLayerPotentialSource.density_discr`
-    of the layer potential source identified by :attr:`where`.
-    """
+        if self.granularity == GRANULARITY_CENTER:
+            name.append("center")
+        elif self.granularity == GRANULARITY_ELEMENT:
+            name.append("panel")
 
+        return "/".join(name)
 
-class QBXSourceStage2(_QBXSource):
-    """A symbolic 'where' specifier for the
-    :attr:`pytential.qbx.QBXLayerPotentialSource.stage2_density_discr`
-    of the layer potential source identified by :attr:`where`.
-    """
 
+def as_dofdesc(desc):
+    if isinstance(desc, DOFDescriptor):
+        return desc
 
-class QBXSourceQuadStage2(_QBXSource):
-    """A symbolic 'where' specifier for the
-    :attr:`pytential.qbx.QBXLayerPotentialSource.quad_stage2_density_discr`
-    of the layer potential source identified by :attr:`where`.
-    """
+    if desc == QBX_SOURCE_STAGE1 \
+            or desc == QBX_SOURCE_STAGE2 \
+            or desc == QBX_SOURCE_QUAD_STAGE2:
+        return DOFDescriptor(discr_stage=desc)
+
+    if desc == GRANULARITY_NODE \
+            or desc == GRANULARITY_CENTER \
+            or desc == GRANULARITY_ELEMENT:
+        return DOFDescriptor(granularity=desc)
+
+    return DOFDescriptor(geometry=desc)
 
 # }}}
 
 
-class cse_scope(cse_scope_base):  # noqa
+class cse_scope(cse_scope_base):  # noqa: N801
     DISCRETIZATION = "pytential_discretization"
 
 
@@ -287,8 +443,9 @@ def make_sym_mv(name, num_components):
     return MultiVector(make_sym_vector(name, num_components))
 
 
-def make_sym_surface_mv(name, ambient_dim, dim, where=None):
-    par_grad = parametrization_derivative_matrix(ambient_dim, dim, where)
+@_deprecate_kwargs('where', 'dofdesc')
+def make_sym_surface_mv(name, ambient_dim, dim, dofdesc=None):
+    par_grad = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     return sum(
             var("%s%d" % (name, i))
@@ -353,17 +510,24 @@ class DiscretizationProperty(Expression):
     further arguments.
     """
 
-    init_arg_names = ("where",)
+    init_arg_names = ("dofdesc",)
 
-    def __init__(self, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, dofdesc=None):
         """
-        :arg where: |where-blurb|
+        :arg dofdesc: |dofdesc-blurb|
         """
 
-        self.where = where
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
-        return (self.where,)
+        return (self.dofdesc,)
 
 
 # {{{ discretization properties
@@ -376,29 +540,31 @@ class QWeight(DiscretizationProperty):
 
 class NodeCoordinateComponent(DiscretizationProperty):
 
-    init_arg_names = ("ambient_axis", "where")
+    init_arg_names = ("ambient_axis", "dofdesc")
 
-    def __init__(self, ambient_axis, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, ambient_axis, dofdesc=None):
         """
-        :arg where: |where-blurb|
+        :arg dofdesc: |dofdesc-blurb|
         """
         self.ambient_axis = ambient_axis
-        DiscretizationProperty.__init__(self, where)
+        DiscretizationProperty.__init__(self, dofdesc)
 
     def __getinitargs__(self):
-        return (self.ambient_axis, self.where)
+        return (self.ambient_axis, self.dofdesc)
 
     mapper_method = intern("map_node_coordinate_component")
 
 
-def nodes(ambient_dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def nodes(ambient_dim, dofdesc=None):
     """Return a :class:`pymbolic.geometric_algebra.MultiVector` of node
     locations.
     """
 
     return MultiVector(
             make_obj_array([
-                NodeCoordinateComponent(i, where)
+                NodeCoordinateComponent(i, dofdesc)
                 for i in range(ambient_dim)]))
 
 
@@ -408,22 +574,24 @@ class NumReferenceDerivative(DiscretizationProperty):
     reference coordinates.
     """
 
-    init_arg_names = ("ref_axes", "operand", "where")
+    init_arg_names = ("ref_axes", "operand", "dofdesc")
 
-    def __new__(cls, ref_axes=None, operand=None, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __new__(cls, ref_axes=None, operand=None, dofdesc=None):
         # If the constructor is handed a multivector object, return an
         # object array of the operator applied to each of the
         # coefficients in the multivector.
 
         if isinstance(operand, np.ndarray):
             def make_op(operand_i):
-                return cls(ref_axes, operand_i, where=where)
+                return cls(ref_axes, operand_i, dofdesc=dofdesc)
 
             return componentwise(make_op, operand)
         else:
             return DiscretizationProperty.__new__(cls)
 
-    def __init__(self, ref_axes, operand, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, ref_axes, operand, dofdesc=None):
         """
         :arg ref_axes: a :class:`tuple` of tuples indicating indices of
             coordinate axes of the reference element to the number of derivatives
@@ -433,7 +601,7 @@ class NumReferenceDerivative(DiscretizationProperty):
 
             May also be a singile integer *i*, which is viewed as equivalent
             to ``((i, 1),)``.
-        :arg where: |where-blurb|
+        :arg dofdesc: |dofdesc-blurb|
         """
 
         if isinstance(ref_axes, int):
@@ -451,15 +619,16 @@ class NumReferenceDerivative(DiscretizationProperty):
         self.ref_axes = ref_axes
 
         self.operand = operand
-        DiscretizationProperty.__init__(self, where)
+        DiscretizationProperty.__init__(self, dofdesc)
 
     def __getinitargs__(self):
-        return (self.ref_axes, self.operand, self.where)
+        return (self.ref_axes, self.operand, self.dofdesc)
 
     mapper_method = intern("map_num_reference_derivative")
 
 
-def reference_jacobian(func, output_dim, dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def reference_jacobian(func, output_dim, dim, dofdesc=None):
     """Return a :class:`np.array` representing the Jacobian of a vector function
     with respect to the reference coordinates.
     """
@@ -468,35 +637,38 @@ def reference_jacobian(func, output_dim, dim, where=None):
     for i in range(output_dim):
         func_component = func[i]
         for j in range(dim):
-            jac[i, j] = NumReferenceDerivative(j, func_component, where)
+            jac[i, j] = NumReferenceDerivative(j, func_component, dofdesc)
 
     return jac
 
 
-def parametrization_derivative_matrix(ambient_dim, dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def parametrization_derivative_matrix(ambient_dim, dim, dofdesc=None):
     """Return a :class:`np.array` representing the derivative of the
     reference-to-global parametrization.
     """
 
     return cse(
             reference_jacobian(
-                [NodeCoordinateComponent(i, where) for i in range(ambient_dim)],
-                ambient_dim, dim, where=where),
+                [NodeCoordinateComponent(i, dofdesc) for i in range(ambient_dim)],
+                ambient_dim, dim, dofdesc=dofdesc),
             "pd_matrix", cse_scope.DISCRETIZATION)
 
 
-def parametrization_derivative(ambient_dim, dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def parametrization_derivative(ambient_dim, dim, dofdesc=None):
     """Return a :class:`pymbolic.geometric_algebra.MultiVector` representing
     the derivative of the reference-to-global parametrization.
     """
 
-    par_grad = parametrization_derivative_matrix(ambient_dim, dim, where)
+    par_grad = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     from pytools import product
     return product(MultiVector(vec) for vec in par_grad.T)
 
 
-def pseudoscalar(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def pseudoscalar(ambient_dim, dim=None, dofdesc=None):
     """
     Same as the outer product of all parametrization derivative columns.
     """
@@ -504,34 +676,37 @@ def pseudoscalar(ambient_dim, dim=None, where=None):
         dim = ambient_dim - 1
 
     return cse(
-            parametrization_derivative(ambient_dim, dim, where)
+            parametrization_derivative(ambient_dim, dim, dofdesc)
             .project_max_grade(),
             "pseudoscalar", cse_scope.DISCRETIZATION)
 
 
-def area_element(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def area_element(ambient_dim, dim=None, dofdesc=None):
     return cse(
-            sqrt(pseudoscalar(ambient_dim, dim, where).norm_squared()),
+            sqrt(pseudoscalar(ambient_dim, dim, dofdesc).norm_squared()),
             "area_element", cse_scope.DISCRETIZATION)
 
 
-def sqrt_jac_q_weight(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def sqrt_jac_q_weight(ambient_dim, dim=None, dofdesc=None):
     return cse(
             sqrt(
-                area_element(ambient_dim, dim, where)
-                * QWeight(where)),
+                area_element(ambient_dim, dim, dofdesc)
+                * QWeight(dofdesc)),
             "sqrt_jac_q_weight", cse_scope.DISCRETIZATION)
 
 
-def normal(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def normal(ambient_dim, dim=None, dofdesc=None):
     """Exterior unit normals."""
 
     # Don't be tempted to add a sign here. As it is, it produces
     # exterior normals for positively oriented curves and surfaces.
 
     pder = (
-            pseudoscalar(ambient_dim, dim, where)
-            / area_element(ambient_dim, dim, where))
+            pseudoscalar(ambient_dim, dim, dofdesc)
+            / area_element(ambient_dim, dim, dofdesc))
     return cse(
             # Dorst Section 3.7.2
             pder << pder.I.inv(),
@@ -539,7 +714,8 @@ def normal(ambient_dim, dim=None, where=None):
             scope=cse_scope.DISCRETIZATION)
 
 
-def mean_curvature(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def mean_curvature(ambient_dim, dim=None, dofdesc=None):
     """(Numerical) mean curvature."""
 
     if dim is None:
@@ -547,16 +723,16 @@ def mean_curvature(ambient_dim, dim=None, where=None):
 
     if ambient_dim == 2 and dim == 1:
         # https://en.wikipedia.org/wiki/Curvature#Local_expressions
-        xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where)
+        xp, yp = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
         xpp, ypp = cse(
-                reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where),
+                reference_jacobian([xp[0], yp[0]], ambient_dim, dim, dofdesc),
                 "p2d_matrix", cse_scope.DISCRETIZATION)
 
         kappa = (xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2)
     elif ambient_dim == 3 and dim == 2:
         # https://en.wikipedia.org/wiki/Mean_curvature#Surfaces_in_3D_space
-        s_op = shape_operator(ambient_dim, dim=dim, where=where)
+        s_op = shape_operator(ambient_dim, dim=dim, dofdesc=dofdesc)
         kappa = -0.5 * sum(s_op[i, i] for i in range(s_op.shape[0]))
     else:
         raise NotImplementedError('not available in {}D for {}D surfaces'
@@ -565,21 +741,23 @@ def mean_curvature(ambient_dim, dim=None, where=None):
     return kappa
 
 
-def first_fundamental_form(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def first_fundamental_form(ambient_dim, dim=None, dofdesc=None):
     if dim is None:
         dim = ambient_dim - 1
 
     if ambient_dim != 3 and dim != 2:
         raise NotImplementedError("only available for surfaces in 3D")
 
-    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
+    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     return cse(
             np.dot(pd_mat.T, pd_mat),
             "fundform1")
 
 
-def second_fundamental_form(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def second_fundamental_form(ambient_dim, dim=None, dofdesc=None):
     """Compute the second fundamental form of a surface. This is in reference
     to the reference-to-global mapping in use for each element.
 
@@ -594,17 +772,17 @@ def second_fundamental_form(ambient_dim, dim=None, where=None):
     if ambient_dim != 3 and dim != 2:
         raise NotImplementedError("only available for surfaces in 3D")
 
-    r = nodes(ambient_dim, where=where).as_vector()
+    r = nodes(ambient_dim, dofdesc=dofdesc).as_vector()
 
     # https://en.wikipedia.org/w/index.php?title=Second_fundamental_form&oldid=821047433#Classical_notation
 
     from functools import partial
-    d = partial(NumReferenceDerivative, where=where)
+    d = partial(NumReferenceDerivative, dofdesc=dofdesc)
     ruu = d(((0, 2),), r)
     ruv = d(((0, 1), (1, 1)), r)
     rvv = d(((1, 2),), r)
 
-    nrml = normal(ambient_dim, dim, where).as_vector()
+    nrml = normal(ambient_dim, dim, dofdesc).as_vector()
 
     ff2_l = cse(np.dot(ruu, nrml), "fundform2_L")
     ff2_m = cse(np.dot(ruv, nrml), "fundform2_M")
@@ -618,7 +796,8 @@ def second_fundamental_form(ambient_dim, dim=None, where=None):
     return result
 
 
-def shape_operator(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def shape_operator(ambient_dim, dim=None, dofdesc=None):
     if dim is None:
         dim = ambient_dim - 1
 
@@ -626,8 +805,8 @@ def shape_operator(ambient_dim, dim=None, where=None):
         raise NotImplementedError("only available for surfaces in 3D")
 
     # https://en.wikipedia.org/w/index.php?title=Differential_geometry_of_surfaces&oldid=833587563
-    (E, F), (F, G) = first_fundamental_form(ambient_dim, dim, where)
-    (e, f), (f, g) = second_fundamental_form(ambient_dim, dim, where)
+    (E, F), (F, G) = first_fundamental_form(ambient_dim, dim, dofdesc)
+    (e, f), (f, g) = second_fundamental_form(ambient_dim, dim, dofdesc)
 
     result = np.zeros((2, 2), dtype=object)
     result[0, 0] = e*G-f*F
@@ -640,7 +819,8 @@ def shape_operator(ambient_dim, dim=None, where=None):
             "shape_operator")
 
 
-def _panel_size(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _panel_size(ambient_dim, dim=None, dofdesc=None):
     # A broken quasi-1D approximation of 1D element size. Do not use.
 
     if dim is None:
@@ -687,12 +867,13 @@ def _small_mat_eigenvalues(mat):
                 "eigenvalue formula for %dx%d matrices" % (m, n))
 
 
+@_deprecate_kwargs('where', 'dofdesc')
 def _equilateral_parametrization_derivative_matrix(ambient_dim, dim=None,
-        where=None):
+        dofdesc=None):
     if dim is None:
         dim = ambient_dim - 1
 
-    pder_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
+    pder_mat = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     # The above procedure works well only when the 'reference' end of the
     # mapping is in equilateral coordinates.
@@ -705,13 +886,15 @@ def _equilateral_parametrization_derivative_matrix(ambient_dim, dim=None,
             "equilateral_pder_mat")
 
 
-def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None,
+@_deprecate_kwargs('where', 'dofdesc')
+def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, dofdesc=None,
         with_elementwise_max=True):
     """Return the largest factor by which the reference-to-global
     mapping stretches the bi-unit (i.e. :math:`[-1,1]`) reference
     element along any axis.
 
-    Returns a DOF vector that is elementwise constant.
+    If *map_elementwise_max* is True, returns a DOF vector that is elementwise
+    constant.
     """
 
     if dim is None:
@@ -727,7 +910,7 @@ def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None,
     # reflects available quadrature resolution in that direction.
 
     equi_pder_mat = _equilateral_parametrization_derivative_matrix(
-            ambient_dim, dim, where)
+            ambient_dim, dim, dofdesc)
 
     # Compute eigenvalues of J^T to compute SVD.
     equi_pder_mat_jtj = cse(
@@ -747,21 +930,22 @@ def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=None,
     result = Max(tuple(stretch_factors))
 
     if with_elementwise_max:
-        result = ElementwiseMax(result, where=where)
+        result = ElementwiseMax(result, dofdesc=dofdesc)
 
     return cse(result, "mapping_max_stretch", cse_scope.DISCRETIZATION)
 
 
-def _max_curvature(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _max_curvature(ambient_dim, dim=None, dofdesc=None):
     # An attempt at a 'max curvature' criterion.
 
     if dim is None:
         dim = ambient_dim - 1
 
     if ambient_dim == 2:
-        return abs(mean_curvature(ambient_dim, dim, where=where))
+        return abs(mean_curvature(ambient_dim, dim, dofdesc=dofdesc))
     elif ambient_dim == 3:
-        shape_op = shape_operator(ambient_dim, dim, where=where)
+        shape_op = shape_operator(ambient_dim, dim, dofdesc=dofdesc)
 
         abs_principal_curvatures = [
                 abs(x) for x in _small_mat_eigenvalues(shape_op)]
@@ -772,22 +956,170 @@ def _max_curvature(ambient_dim, dim=None, where=None):
                 "dimensions" % ambient_dim)
 
 
-def _scaled_max_curvature(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def _scaled_max_curvature(ambient_dim, dim=None, dofdesc=None):
     """An attempt at a unit-less, scale-invariant quantity that characterizes
     'how much curviness there is on an element'. Values seem to hover around 1
     on typical meshes. Empirical evidence suggests that elements exceeding
     a threshold of about 0.8-1 will have high QBX truncation error.
     """
 
-    return _max_curvature(ambient_dim, dim, where=where) * \
-            _simplex_mapping_max_stretch_factor(ambient_dim, dim, where=where,
+    return _max_curvature(ambient_dim, dim, dofdesc=dofdesc) * \
+            _simplex_mapping_max_stretch_factor(ambient_dim, dim,
+                    dofdesc=dofdesc,
                     with_elementwise_max=False)
 
 # }}}
 
 
+# {{{ qbx-specific geometry
+
+def _expansion_radii_factor(ambient_dim, dim):
+    if dim is None:
+        dim = ambient_dim - 1
+
+    dim_fudge_factor = 0.5 if dim == 2 else 1.0
+    return 0.5 * dim_fudge_factor
+
+
+@_deprecate_kwargs('where', 'dofdesc')
+def _quad_resolution(ambient_dim, dim=None, granularity=None, dofdesc=None):
+    """This measures the quadrature resolution across the
+    mesh. In a 1D uniform mesh of uniform 'parametrization speed', it
+    should be the same as the panel length.
+
+    In multiple dimensions (i.e. with multiple quadrature resolutions
+    depending on direction), this measure returns the coarsest of these resolution,
+    is invariant with respect to rotation of the global coordinate
+    system, and invariant with respect to vertex ordering of the reference
+    element.
+    """
+
+    from_dd = as_dofdesc(dofdesc)
+    to_dd = from_dd.copy(granularity=granularity)
+
+    stretch = _simplex_mapping_max_stretch_factor(ambient_dim,
+            dim=dim, dofdesc=from_dd)
+    return interp(from_dd, to_dd, stretch)
+
+
+@_deprecate_kwargs('where', 'dofdesc')
+def _source_danger_zone_radii(ambient_dim, dim=None,
+        granularity=None, dofdesc=None):
+    dofdesc = as_dofdesc(dofdesc)
+    if dofdesc.discr_stage is None:
+        dofdesc = dofdesc.copy(discr_stage=QBX_SOURCE_STAGE2)
+
+    # This should be the expression of the expansion radii, but
+    #
+    # - in reference to the stage 2 discretization
+    # - mutliplied by 0.75 because
+    #
+    #   - Setting this equal to the expansion radii ensures that *every*
+    #     stage 2 element will be refined, which is wasteful.
+    #     (so this needs to be smaller than that)
+    #   - Setting this equal to half the expansion radius will not provide
+    #     a refinement 'buffer layer' at a 2x coarsening fringe.
+
+    factor = 0.75 * _expansion_radii_factor(ambient_dim, dim)
+    return factor * _quad_resolution(ambient_dim, dim=dim,
+            granularity=granularity, dofdesc=dofdesc)
+
+
+@_deprecate_kwargs('where', 'dofdesc')
+def _close_target_tunnel_radii(ambient_dim, dim=None,
+        granularity=None, dofdesc=None):
+    factor = 0.5 * _expansion_radii_factor(ambient_dim, dim)
+
+    return factor * _quad_resolution(ambient_dim, dim=dim,
+            granularity=granularity, dofdesc=dofdesc)
+
+
+@_deprecate_kwargs('where', 'dofdesc')
+def expansion_radii(ambient_dim, dim=None, granularity=None, dofdesc=None):
+    factor = _expansion_radii_factor(ambient_dim, dim)
+
+    return cse(factor * _quad_resolution(ambient_dim, dim=dim,
+        granularity=granularity, dofdesc=dofdesc),
+        "expansion_radii",
+        cse_scope.DISCRETIZATION)
+
+
+@_deprecate_kwargs('where', 'dofdesc')
+def expansion_centers(ambient_dim, side, dim=None, dofdesc=None):
+    x = nodes(ambient_dim, dofdesc=dofdesc)
+    normals = normal(ambient_dim, dim=dim, dofdesc=dofdesc)
+    radii = expansion_radii(ambient_dim, dim=dim,
+            granularity=GRANULARITY_NODE, dofdesc=dofdesc)
+
+    centers = x + side * radii * normals
+    return cse(centers.as_vector(),
+            "expansion_centers",
+            cse_scope.DISCRETIZATION)
+
+
+@_deprecate_kwargs('where', 'dofdesc')
+def h_max(ambient_dim, dim=None, dofdesc=None):
+    """Defines a maximum element size in the discretization."""
+
+    dofdesc = as_dofdesc(dofdesc).copy(granularity=GRANULARITY_ELEMENT)
+    r = _quad_resolution(ambient_dim, dim=dim, dofdesc=dofdesc)
+
+    return cse(NodeMax(r),
+            "h_max",
+            cse_scope.DISCRETIZATION)
+
+
+@_deprecate_kwargs('where', 'dofdesc')
+def weights_and_area_elements(ambient_dim, dim=None, dofdesc=None):
+    """Combines :func:`area_element` and :class:`QWeight`."""
+
+    area = area_element(ambient_dim, dim=dim, dofdesc=dofdesc)
+    return cse(area * QWeight(dofdesc=dofdesc),
+            "weights_area_elements",
+            cse_scope.DISCRETIZATION)
+
+# }}}
+
+
 # {{{ operators
 
+class Interpolation(Expression):
+    """Interpolate quantity from a DOF described by *from_dd* to a DOF
+    described by *to_dd*."""
+
+    init_arg_names = ("from_dd", "to_dd", "operand")
+
+    def __new__(cls, from_dd, to_dd, operand):
+        from_dd = as_dofdesc(from_dd)
+        to_dd = as_dofdesc(to_dd)
+
+        if from_dd == to_dd:
+            return operand
+
+        if isinstance(operand, np.ndarray):
+            def make_op(operand_i):
+                return cls(from_dd, to_dd, operand_i)
+
+            return componentwise(make_op, operand)
+        else:
+            return Expression.__new__(cls)
+
+    def __init__(self, from_dd, to_dd, operand):
+        self.from_dd = as_dofdesc(from_dd)
+        self.to_dd = as_dofdesc(to_dd)
+        self.operand = operand
+
+    def __getinitargs__(self):
+        return (self.from_dd, self.to_dd, self.operand)
+
+    mapper_method = intern("map_interpolation")
+
+
+def interp(from_dd, to_dd, operand):
+    return Interpolation(from_dd, to_dd, operand)
+
+
 class SingleScalarOperandExpression(Expression):
 
     init_arg_names = ("operand",)
@@ -824,38 +1156,47 @@ class NodeMax(SingleScalarOperandExpression):
     mapper_method = "map_node_max"
 
 
-def integral(ambient_dim, dim, operand, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def integral(ambient_dim, dim, operand, dofdesc=None):
     """A volume integral of *operand*."""
 
     return NodeSum(
-            area_element(ambient_dim, dim, where)
-            * QWeight(where)
+            area_element(ambient_dim, dim, dofdesc)
+            * QWeight(dofdesc)
             * operand)
 
 
 class SingleScalarOperandExpressionWithWhere(Expression):
 
-    init_arg_names = ("operand", "where")
+    init_arg_names = ("operand", "dofdesc")
 
-    def __new__(cls, operand=None, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __new__(cls, operand=None, dofdesc=None):
         # If the constructor is handed a multivector object, return an
         # object array of the operator applied to each of the
         # coefficients in the multivector.
 
         if isinstance(operand, (np.ndarray, MultiVector)):
             def make_op(operand_i):
-                return cls(operand_i, where)
+                return cls(operand_i, dofdesc)
 
             return componentwise(make_op, operand)
         else:
             return Expression.__new__(cls)
 
-    def __init__(self, operand, where=None):
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, operand, dofdesc=None):
         self.operand = operand
-        self.where = where
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
-        return (self.operand, self.where)
+        return (self.operand, self.dofdesc)
 
 
 class ElementwiseSum(SingleScalarOperandExpressionWithWhere):
@@ -887,54 +1228,71 @@ class Ones(Expression):
     discretization.
     """
 
-    init_arg_names = ("where",)
+    init_arg_names = ("dofdesc",)
 
-    def __init__(self, where=None):
-        self.where = where
+    @_deprecate_kwargs('where', 'dofdesc')
+    def __init__(self, dofdesc=None):
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
-        return (self.where,)
+        return (self.dofdesc,)
 
     mapper_method = intern("map_ones")
 
 
-def ones_vec(dim, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def ones_vec(dim, dofdesc=None):
     from pytools.obj_array import make_obj_array
     return MultiVector(
-                make_obj_array(dim*[Ones(where)]))
+                make_obj_array(dim*[Ones(dofdesc)]))
 
 
-def area(ambient_dim, dim, where=None):
-    return cse(integral(ambient_dim, dim, Ones(where), where), "area",
+@_deprecate_kwargs('where', 'dofdesc')
+def area(ambient_dim, dim, dofdesc=None):
+    return cse(integral(ambient_dim, dim, Ones(dofdesc), dofdesc), "area",
             cse_scope.DISCRETIZATION)
 
 
-def mean(ambient_dim, dim, operand, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def mean(ambient_dim, dim, operand, dofdesc=None):
     return (
-            integral(ambient_dim, dim, operand, where)
-            / area(ambient_dim, dim, where))
+            integral(ambient_dim, dim, operand, dofdesc)
+            / area(ambient_dim, dim, dofdesc))
 
 
 class IterativeInverse(Expression):
 
-    init_arg_names = ("expression", "rhs", "variable_name", "extra_vars", "where")
+    init_arg_names = ("expression", "rhs", "variable_name", "extra_vars", "dofdesc")
 
+    @_deprecate_kwargs('where', 'dofdesc')
     def __init__(self, expression, rhs, variable_name, extra_vars={},
-            where=None):
+            dofdesc=None):
         self.expression = expression
         self.rhs = rhs
         self.variable_name = variable_name
         self.extra_vars = extra_vars
-        self.where = where
+        self.dofdesc = as_dofdesc(dofdesc)
+
+    @property
+    def where(self):
+        warn('`where` is deprecated. use `dofdesc` instead.',
+             DeprecationWarning, stacklevel=2)
+        return self.dofdesc
 
     def __getinitargs__(self):
         return (self.expression, self.rhs, self.variable_name,
-                self.extra_vars, self.where)
+                self.extra_vars, self.dofdesc)
 
     def get_hash(self):
         return hash((self.__class__,) + (self.expression,
             self.rhs, self.variable_name,
-            frozenset(six.iteritems(self.extra_vars)), self.where))
+            frozenset(six.iteritems(self.extra_vars)), self.dofdesc))
 
     mapper_method = intern("map_inverse")
 
@@ -1022,31 +1380,6 @@ class _NoArgSentinel(object):
     pass
 
 
-class Interpolation(Expression):
-    """Interpolate quantity from *source* to *target* discretization."""
-
-    init_arg_names = ("source", "target", "operand")
-
-    def __new__(cls, source, target, operand):
-        if isinstance(operand, np.ndarray):
-            def make_op(operand_i):
-                return cls(source, target, operand_i)
-
-            return componentwise(make_op, operand)
-        else:
-            return Expression.__new__(cls)
-
-    def __init__(self, source, target, operand):
-        self.source = source
-        self.target = target
-        self.operand = operand
-
-    def __getinitargs__(self):
-        return (self.source, self.target, self.operand)
-
-    mapper_method = intern("map_interpolation")
-
-
 class IntG(Expression):
     r"""
     .. math::
@@ -1169,8 +1502,8 @@ class IntG(Expression):
         self.kernel = kernel
         self.density = density
         self.qbx_forced_limit = qbx_forced_limit
-        self.source = source
-        self.target = target
+        self.source = as_dofdesc(source)
+        self.target = as_dofdesc(target)
         self.kernel_arguments = kernel_arguments
 
     def copy(self, kernel=None, density=None, qbx_forced_limit=_NoArgSentinel,
@@ -1179,8 +1512,8 @@ class IntG(Expression):
         density = density or self.density
         if qbx_forced_limit is _NoArgSentinel:
             qbx_forced_limit = self.qbx_forced_limit
-        source = source or self.source
-        target = target or self.target
+        source = as_dofdesc(source or self.source)
+        target = as_dofdesc(target or self.target)
         kernel_arguments = kernel_arguments or self.kernel_arguments
         return type(self)(kernel, density, qbx_forced_limit, source, target,
                 kernel_arguments)
@@ -1299,11 +1632,11 @@ def int_g_dsource(ambient_dim, dsource, kernel, density,
 # {{{ geometric calculus
 
 
-class _unspecified:  # noqa
+class _unspecified:  # noqa: N801
     pass
 
 
-def S(kernel, density,  # noqa
+def S(kernel, density,
         qbx_forced_limit=_unspecified, source=None, target=None,
         kernel_arguments=None, **kwargs):
 
@@ -1316,10 +1649,11 @@ def S(kernel, density,  # noqa
             kernel_arguments, **kwargs)
 
 
-def tangential_derivative(ambient_dim, operand, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def tangential_derivative(ambient_dim, operand, dim=None, dofdesc=None):
     pder = (
-            pseudoscalar(ambient_dim, dim, where)
-            / area_element(ambient_dim, dim, where))
+            pseudoscalar(ambient_dim, dim, dofdesc)
+            / area_element(ambient_dim, dim, dofdesc))
 
     # FIXME: Should be formula (3.25) in Dorst et al.
     d = Derivative()
@@ -1327,15 +1661,16 @@ def tangential_derivative(ambient_dim, operand, dim=None, where=None):
             (d.dnabla(ambient_dim) * d(operand)) >> pder)
 
 
-def normal_derivative(ambient_dim, operand, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def normal_derivative(ambient_dim, operand, dim=None, dofdesc=None):
     d = Derivative()
     return d.resolve(
-            (normal(ambient_dim, dim, where).scalar_product(d.dnabla(ambient_dim)))
+            (normal(ambient_dim, dim, dofdesc).scalar_product(d.dnabla(ambient_dim)))
             * d(operand))
 
 
-def Sp(kernel, *args, **kwargs):  # noqa
-    where = kwargs.get("target")
+def Sp(kernel, *args, **kwargs):
+    dofdesc = kwargs.get("target")
     if "qbx_forced_limit" not in kwargs:
         warn("not specifying qbx_forced_limit on call to 'Sp' is deprecated, "
                 "defaulting to 'avg'", DeprecationWarning, stacklevel=2)
@@ -1353,10 +1688,10 @@ def Sp(kernel, *args, **kwargs):  # noqa
     return normal_derivative(
             ambient_dim,
             S(kernel, *args, **kwargs),
-            dim=dim, where=where)
+            dim=dim, dofdesc=dofdesc)
 
 
-def Spp(kernel, *args, **kwargs):  # noqa
+def Spp(kernel, *args, **kwargs):
     ambient_dim = kwargs.get("ambient_dim")
     from sumpy.kernel import Kernel
     if ambient_dim is None and isinstance(kernel, Kernel):
@@ -1366,14 +1701,14 @@ def Spp(kernel, *args, **kwargs):  # noqa
                 "the kernel, or directly")
     dim = kwargs.pop("dim", None)
 
-    where = kwargs.get("target")
+    dofdesc = kwargs.get("target")
     return normal_derivative(
             ambient_dim,
             Sp(kernel, *args, **kwargs),
-            dim=dim, where=where)
+            dim=dim, dofdesc=dofdesc)
 
 
-def D(kernel, *args, **kwargs):  # noqa
+def D(kernel, *args, **kwargs):
     ambient_dim = kwargs.get("ambient_dim")
     from sumpy.kernel import Kernel
     if ambient_dim is None and isinstance(kernel, Kernel):
@@ -1383,7 +1718,7 @@ def D(kernel, *args, **kwargs):  # noqa
                 "the kernel, or directly")
     dim = kwargs.pop("dim", None)
 
-    where = kwargs.get("source")
+    dofdesc = kwargs.get("source")
 
     if "qbx_forced_limit" not in kwargs:
         warn("not specifying qbx_forced_limit on call to 'D' is deprecated, "
@@ -1392,11 +1727,11 @@ def D(kernel, *args, **kwargs):  # noqa
 
     return int_g_dsource(
             ambient_dim,
-            normal(ambient_dim, dim, where),
+            normal(ambient_dim, dim, dofdesc),
             kernel, *args, **kwargs).xproject(0)
 
 
-def Dp(kernel, *args, **kwargs):  # noqa
+def Dp(kernel, *args, **kwargs):
     ambient_dim = kwargs.get("ambient_dim")
     from sumpy.kernel import Kernel
     if ambient_dim is None and isinstance(kernel, Kernel):
@@ -1413,7 +1748,7 @@ def Dp(kernel, *args, **kwargs):  # noqa
     return normal_derivative(
             ambient_dim,
             D(kernel, *args, **kwargs),
-            dim=dim, where=target)
+            dim=dim, dofdesc=target)
 
 # }}}
 
@@ -1424,15 +1759,16 @@ def Dp(kernel, *args, **kwargs):  # noqa
 
 # {{{ conventional vector calculus
 
-def tangential_onb(ambient_dim, dim=None, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def tangential_onb(ambient_dim, dim=None, dofdesc=None):
     """Return a matrix of shape ``(ambient_dim, dim)`` with orthogonal columns
-    spanning the tangential space of the surface of *where*.
+    spanning the tangential space of the surface of *dofdesc*.
     """
 
     if dim is None:
         dim = ambient_dim - 1
 
-    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
+    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, dofdesc)
 
     # {{{ Gram-Schmidt
 
@@ -1451,30 +1787,34 @@ def tangential_onb(ambient_dim, dim=None, where=None):
     return orth_pd_mat
 
 
-def xyz_to_tangential(xyz_vec, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def xyz_to_tangential(xyz_vec, dofdesc=None):
     ambient_dim = len(xyz_vec)
-    tonb = tangential_onb(ambient_dim, where=where)
+    tonb = tangential_onb(ambient_dim, dofdesc=dofdesc)
     return make_obj_array([
         np.dot(tonb[:, i], xyz_vec)
         for i in range(ambient_dim - 1)
         ])
 
 
-def tangential_to_xyz(tangential_vec, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def tangential_to_xyz(tangential_vec, dofdesc=None):
     ambient_dim = len(tangential_vec) + 1
-    tonb = tangential_onb(ambient_dim, where=where)
+    tonb = tangential_onb(ambient_dim, dofdesc=dofdesc)
     return sum(
         tonb[:, i] * tangential_vec[i]
         for i in range(ambient_dim - 1))
 
 
-def project_to_tangential(xyz_vec, where=None):
+@_deprecate_kwargs('where', 'dofdesc')
+def project_to_tangential(xyz_vec, dofdesc=None):
     return tangential_to_xyz(
-            cse(xyz_to_tangential(xyz_vec, where), where))
+            cse(xyz_to_tangential(xyz_vec, dofdesc), dofdesc))
 
 
-def n_dot(vec, where=None):
-    nrm = normal(len(vec), where).as_vector()
+@_deprecate_kwargs('where', 'dofdesc')
+def n_dot(vec, dofdesc=None):
+    nrm = normal(len(vec), dofdesc).as_vector()
 
     return np.dot(nrm, vec)
 
@@ -1491,8 +1831,9 @@ def cross(vec_a, vec_b):
         for i in range(3)])
 
 
-def n_cross(vec, where=None):
-    return cross(normal(3, where).as_vector(), vec)
+@_deprecate_kwargs('where', 'dofdesc')
+def n_cross(vec, dofdesc=None):
+    return cross(normal(3, dofdesc).as_vector(), vec)
 
 
 def div(vec):
diff --git a/pytential/unregularized.py b/pytential/unregularized.py
index 4bd67ae97fe1206442bfdad357ec8cbc9ef0779f..2ff748d18d31303fb1f73998bdd72851ca66e47e 100644
--- a/pytential/unregularized.py
+++ b/pytential/unregularized.py
@@ -87,37 +87,14 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
             expansion_factory = DefaultExpansionFactory()
         self.expansion_factory = expansion_factory
 
-        self.debug = debug
-
     @memoize_method
     def weights_and_area_elements(self):
-        import pytential.symbolic.primitives as p
-        from pytential.symbolic.execution import bind
+        from pytential import bind, sym
         with cl.CommandQueue(self.cl_context) as queue:
-            # quad_stage2_density_discr is not guaranteed to be usable for
-            # interpolation/differentiation. Use density_discr to find
-            # area element instead, then upsample that.
-
-            area_element = self.resampler(
-                    queue,
-                    bind(
-                        self.density_discr,
-                        p.area_element(self.ambient_dim, self.dim)
-                        )(queue))
-
-            qweight = bind(self.quad_stage2_density_discr, p.QWeight())(queue)
+            waa = bind(self,
+                    sym.weights_and_area_elements(self.ambient_dim))(queue)
 
-            return (area_element.with_queue(queue)*qweight).with_queue(None)
-
-    @property
-    def quad_stage2_density_discr(self):
-        return self.density_discr
-
-    def resampler(self, queue, f):
-        return f
-
-    def with_refinement(self):
-        raise NotImplementedError
+            return waa.with_queue(None)
 
     def copy(
             self,
@@ -182,7 +159,7 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
 
             evt, output_for_each_kernel = p2p(queue,
                     target_discr.nodes(),
-                    self.quad_stage2_density_discr.nodes(),
+                    self.density_discr.nodes(),
                     [strengths], **kernel_args)
 
             result.append((o.name, output_for_each_kernel[o.kernel_index]))
@@ -223,7 +200,6 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
                 self.debug)
 
     def exec_compute_potential_insn_fmm(self, queue, insn, bound_expr, evaluate):
-
         # {{{ gather unique target discretizations used
 
         target_name_to_index = {}
@@ -236,7 +212,7 @@ class UnregularizedLayerPotentialSource(LayerPotentialSourceBase):
                 continue
 
             target_name_to_index[o.target_name] = len(targets)
-            targets.append(bound_expr.places[o.target_name])
+            targets.append(bound_expr.places.get_geometry(o.target_name))
 
         targets = tuple(targets)
 
@@ -365,11 +341,11 @@ class _FMMGeometryData(object):
 
     @property
     def coord_dtype(self):
-        return self.lpot_source.quad_stage2_density_discr.nodes().dtype
+        return self.lpot_source.density_discr.nodes().dtype
 
     @property
     def ambient_dim(self):
-        return self.lpot_source.quad_stage2_density_discr.ambient_dim
+        return self.lpot_source.density_discr.ambient_dim
 
     @memoize_method
     def traversal(self):
@@ -392,7 +368,7 @@ class _FMMGeometryData(object):
         target_info = self.target_info()
 
         with cl.CommandQueue(self.cl_context) as queue:
-            nsources = lpot_src.quad_stage2_density_discr.nnodes
+            nsources = lpot_src.density_discr.nnodes
             nparticles = nsources + target_info.ntargets
 
             refine_weights = cl.array.zeros(queue, nparticles, dtype=np.int32)
@@ -402,7 +378,7 @@ class _FMMGeometryData(object):
             MAX_LEAF_REFINE_WEIGHT = 32  # noqa
 
             tree, _ = code_getter.build_tree(queue,
-                    particles=lpot_src.quad_stage2_density_discr.nodes(),
+                    particles=lpot_src.density_discr.nodes(),
                     targets=target_info.targets,
                     max_leaf_refine_weight=MAX_LEAF_REFINE_WEIGHT,
                     refine_weights=refine_weights,
diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py
index af340805f7af66c381eae1ce64b8fe6fb6f3ca0d..ca276959ab3ab51b6bec494059a97b0952894283 100644
--- a/test/test_global_qbx.py
+++ b/test/test_global_qbx.py
@@ -42,6 +42,7 @@ from meshmode.mesh.generation import (  # noqa
         make_curve_mesh, generate_icosphere, generate_torus)
 from extra_curve_data import horseshoe
 
+from pytential import bind, sym
 
 import logging
 logger = logging.getLogger(__name__)
@@ -114,19 +115,26 @@ def run_source_refinement_test(ctx_factory, mesh, order, helmholtz_k=None):
                 cl_ctx, TreeCodeContainer(cl_ctx)).get_wrangler(queue),
             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.quad_stage2_density_discr.nodes().get(queue)
-    int_centers = get_centers_on_side(lpot_source, -1)
+
+    int_centers = bind(lpot_source,
+        sym.expansion_centers(lpot_source.ambient_dim, -1))(queue)
     int_centers = np.array([axis.get(queue) for axis in int_centers])
-    ext_centers = get_centers_on_side(lpot_source, +1)
+    ext_centers = bind(lpot_source,
+        sym.expansion_centers(lpot_source.ambient_dim, +1))(queue)
     ext_centers = np.array([axis.get(queue) for axis in ext_centers])
-    expansion_radii = lpot_source._expansion_radii("nsources").get(queue)
-    quad_res = lpot_source._coarsest_quad_resolution("npanels").get(queue)
-    source_danger_zone_radii = \
-            lpot_source._source_danger_zone_radii("npanels").get(queue)
+
+    expansion_radii = bind(lpot_source,
+        sym.expansion_radii(lpot_source.ambient_dim))(queue).get()
+    source_danger_zone_radii = bind(lpot_source, sym._source_danger_zone_radii(
+        lpot_source.ambient_dim,
+        dofdesc=sym.GRANULARITY_ELEMENT))(queue).get()
+
+    quad_res = bind(lpot_source, sym._quad_resolution(
+        lpot_source.ambient_dim,
+        dofdesc=sym.GRANULARITY_ELEMENT))(queue)
 
     # {{{ check if satisfying criteria
 
@@ -258,8 +266,8 @@ def test_target_association(ctx_factory, 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)
-    tunnel_radius = \
-            lpot_source._close_target_tunnel_radius("nsources").with_queue(queue)
+    tunnel_radius = bind(lpot_source,
+        sym._close_target_tunnel_radii(lpot_source.ambient_dim))(queue)
 
     def targets_from_sources(sign, dist):
         from pytential import sym, bind
@@ -316,8 +324,9 @@ def test_target_association(ctx_factory, curve_name, curve_f, nelements,
             target_association_tolerance=1e-10)
         .get(queue=queue))
 
-    expansion_radii = lpot_source._expansion_radii("ncenters").get(queue)
-
+    expansion_radii = bind(lpot_source, sym.expansion_radii(
+        lpot_source.ambient_dim,
+        granularity=sym.GRANULARITY_CENTER))(queue).get()
     surf_targets = np.array(
             [axis.get(queue) for axis in lpot_source.density_discr.nodes()])
     int_targets = np.array([axis.get(queue) for axis in int_targets.nodes()])
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index c99f3073c8a9be4eb535153411eca05ad470e6dc..2f2175aeb078fe2eb6a46cbfc8606e6a931da524 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -486,12 +486,13 @@ def test_3d_jump_relations(ctx_factory, relation, visualize=False):
         bound_jump_identity = bind(qbx, jump_identity_sym)
         jump_identity = bound_jump_identity(queue, density=density)
 
+        h_max = bind(qbx, sym.h_max(qbx.ambient_dim))(queue)
         err = (
                 norm(qbx, queue, jump_identity, np.inf)
                 / norm(qbx, queue, density, np.inf))
-        print("ERROR", qbx.h_max, err)
+        print("ERROR", h_max, err)
 
-        eoc_rec.add_data_point(qbx.h_max, err)
+        eoc_rec.add_data_point(h_max, err)
 
         # {{{ visualization
 
diff --git a/test/test_layer_pot_eigenvalues.py b/test/test_layer_pot_eigenvalues.py
index fdfdc93738d823dae9bcf24828362b54b4bc84ff..5737009ca3d30daa14dd56c39fd5446ea0dde795 100644
--- a/test/test_layer_pot_eigenvalues.py
+++ b/test/test_layer_pot_eigenvalues.py
@@ -119,8 +119,8 @@ def test_ellipse_eigenvalues(ctx_factory, ellipse_aspect, mode_nr, qbx_order,
         if 0:
             # plot geometry, centers, normals
 
-            from pytential.qbx.utils import get_centers_on_side
-            centers = get_centers_on_side(qbx, 1)
+            centers = bind(qbx,
+                    sym.expansion_centers(qbx.ambient_dim, +1))(queue)
 
             nodes_h = nodes.get()
             centers_h = [centers[0].get(), centers[1].get()]
@@ -164,10 +164,11 @@ def test_ellipse_eigenvalues(ctx_factory, ellipse_aspect, mode_nr, qbx_order,
             pt.legend()
             pt.show()
 
+        h_max = bind(qbx, sym.h_max(qbx.ambient_dim))(queue)
         s_err = (
                 norm(density_discr, queue, s_sigma - s_sigma_ref)
                 / norm(density_discr, queue, s_sigma_ref))
-        s_eoc_rec.add_data_point(qbx.h_max, s_err)
+        s_eoc_rec.add_data_point(h_max, s_err)
 
         # }}}
 
@@ -198,7 +199,7 @@ def test_ellipse_eigenvalues(ctx_factory, ellipse_aspect, mode_nr, qbx_order,
         d_err = (
                 norm(density_discr, queue, d_sigma - d_sigma_ref)
                 / d_ref_norm)
-        d_eoc_rec.add_data_point(qbx.h_max, d_err)
+        d_eoc_rec.add_data_point(h_max, d_err)
 
         # }}}
 
@@ -217,7 +218,7 @@ def test_ellipse_eigenvalues(ctx_factory, ellipse_aspect, mode_nr, qbx_order,
             sp_err = (
                     norm(density_discr, queue, sp_sigma - sp_sigma_ref)
                     / norm(density_discr, queue, sigma))
-            sp_eoc_rec.add_data_point(qbx.h_max, sp_err)
+            sp_eoc_rec.add_data_point(h_max, sp_err)
 
             # }}}
 
@@ -313,7 +314,9 @@ def test_sphere_eigenvalues(ctx_factory, mode_m, mode_n, qbx_order,
         s_sigma_op = bind(qbx, sym.S(lap_knl, sym.var("sigma"), qbx_forced_limit=+1))
         s_sigma = s_sigma_op(queue=queue, sigma=ymn)
         s_eigval = 1/(2*mode_n + 1)
-        s_eoc_rec.add_data_point(qbx.h_max, rel_err(s_sigma, s_eigval*ymn))
+
+        h_max = bind(qbx, sym.h_max(qbx.ambient_dim))(queue)
+        s_eoc_rec.add_data_point(h_max, rel_err(s_sigma, s_eigval*ymn))
 
         # }}}
 
@@ -323,7 +326,7 @@ def test_sphere_eigenvalues(ctx_factory, mode_m, mode_n, qbx_order,
                 sym.D(lap_knl, sym.var("sigma"), qbx_forced_limit="avg"))
         d_sigma = d_sigma_op(queue=queue, sigma=ymn)
         d_eigval = -1/(2*(2*mode_n + 1))
-        d_eoc_rec.add_data_point(qbx.h_max, rel_err(d_sigma, d_eigval*ymn))
+        d_eoc_rec.add_data_point(h_max, rel_err(d_sigma, d_eigval*ymn))
 
         # }}}
 
@@ -333,7 +336,8 @@ def test_sphere_eigenvalues(ctx_factory, mode_m, mode_n, qbx_order,
                  sym.Sp(lap_knl, sym.var("sigma"), qbx_forced_limit="avg"))
         sp_sigma = sp_sigma_op(queue=queue, sigma=ymn)
         sp_eigval = -1/(2*(2*mode_n + 1))
-        sp_eoc_rec.add_data_point(qbx.h_max, rel_err(sp_sigma, sp_eigval*ymn))
+
+        sp_eoc_rec.add_data_point(h_max, rel_err(sp_sigma, sp_eigval*ymn))
 
         # }}}
 
@@ -343,7 +347,8 @@ def test_sphere_eigenvalues(ctx_factory, mode_m, mode_n, qbx_order,
                 sym.Dp(lap_knl, sym.var("sigma"), qbx_forced_limit="avg"))
         dp_sigma = dp_sigma_op(queue=queue, sigma=ymn)
         dp_eigval = -(mode_n*(mode_n+1))/(2*mode_n + 1)
-        dp_eoc_rec.add_data_point(qbx.h_max, rel_err(dp_sigma, dp_eigval*ymn))
+
+        dp_eoc_rec.add_data_point(h_max, rel_err(dp_sigma, dp_eigval*ymn))
 
         # }}}
 
diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py
index fa8b7a8f88fe91c647dfb11b8867ec1fb02e7b49..d6d426519e42f1177de0bcf85cf7b8553e8154b5 100644
--- a/test/test_layer_pot_identity.py
+++ b/test/test_layer_pot_identity.py
@@ -382,7 +382,8 @@ def test_identity_convergence(ctx_factory,  case, visualize=False):
         linf_error_norm = norm(density_discr, queue, error, p=np.inf)
         print("--->", key, linf_error_norm)
 
-        eoc_rec.add_data_point(qbx.h_max, linf_error_norm)
+        h_max = bind(qbx, sym.h_max(qbx.ambient_dim))(queue)
+        eoc_rec.add_data_point(h_max, linf_error_norm)
 
         if visualize:
             from meshmode.discretization.visualization import make_visualizer
diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py
index e8a063ca3b94f515e20ee9b130893f0f23243dc3..1707a666026fa639984e96b59daebbd017b61639 100644
--- a/test/test_linalg_proxy.py
+++ b/test/test_linalg_proxy.py
@@ -22,15 +22,13 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
-import os
-import time
-
 import numpy as np
 import numpy.linalg as la
 
 import pyopencl as cl
 from pyopencl.array import to_device
 
+from pytential import bind, sym
 from sumpy.tools import BlockIndexRanges
 from meshmode.mesh.generation import ( # noqa
         ellipse, NArmedStarfish, generate_torus, make_curve_mesh)
@@ -79,32 +77,17 @@ def _build_qbx_discr(queue,
 def _build_block_index(discr,
                        nblks=10,
                        factor=1.0,
-                       method='elements',
                        use_tree=True):
 
-    from pytential.linalg.proxy import (
-            partition_by_nodes, partition_by_elements)
-
-    if method == 'elements':
-        factor = 1.0
+    from pytential.linalg.proxy import partition_by_nodes
 
-    if method == 'nodes':
-        nnodes = discr.nnodes
-    else:
-        nnodes = discr.mesh.nelements
+    nnodes = discr.nnodes
     max_particles_in_box = nnodes // nblks
 
     # create index ranges
-    if method == 'nodes':
-        indices = partition_by_nodes(discr,
-                                     use_tree=use_tree,
-                                     max_nodes_in_box=max_particles_in_box)
-    elif method == 'elements':
-        indices = partition_by_elements(discr,
-                                        use_tree=use_tree,
-                                        max_elements_in_box=max_particles_in_box)
-    else:
-        raise ValueError('unknown method: {}'.format(method))
+    indices = partition_by_nodes(discr,
+                                 use_tree=use_tree,
+                                 max_nodes_in_box=max_particles_in_box)
 
     # randomly pick a subset of points
     if abs(factor - 1.0) > 1.0e-14:
@@ -179,72 +162,23 @@ def _plot_partition_indices(queue, discr, indices, **kwargs):
         vis = make_visualizer(queue, discr, 10)
 
         filename = "test_partition_{0}_{1}_{3}d_{2}.png".format(*args)
-        if os.path.isfile(filename):
-            os.remove(filename)
-
         vis.write_vtk_file(filename, [
             ("marker", cl.array.to_device(queue, marker))
             ])
 
 
-@pytest.mark.parametrize("method", ["nodes", "elements"])
 @pytest.mark.parametrize("use_tree", [True, False])
 @pytest.mark.parametrize("ndim", [2, 3])
-def test_partition_points(ctx_factory, method, use_tree, ndim, visualize=False):
+def test_partition_points(ctx_factory, use_tree, ndim, visualize=False):
     ctx = ctx_factory()
     queue = cl.CommandQueue(ctx)
 
     qbx = _build_qbx_discr(queue, ndim=ndim)
     _build_block_index(qbx.density_discr,
-                       method=method,
                        use_tree=use_tree,
                        factor=0.6)
 
 
-@pytest.mark.parametrize("use_tree", [True, False])
-@pytest.mark.parametrize("ndim", [2, 3])
-def test_partition_coarse(ctx_factory, use_tree, ndim, visualize=False):
-    ctx = ctx_factory()
-    queue = cl.CommandQueue(ctx)
-
-    qbx = _build_qbx_discr(queue, ndim=ndim)
-    srcindices = _build_block_index(qbx.density_discr,
-            method="elements", use_tree=use_tree)
-
-    if visualize:
-        discr = qbx.resampler.from_discr
-        _plot_partition_indices(queue, discr, srcindices,
-                method="elements", use_tree=use_tree, pid="stage1")
-
-    from pytential.linalg.proxy import partition_from_coarse
-    resampler = qbx.direct_resampler
-
-    t_start = time.time()
-    srcindices_ = partition_from_coarse(resampler, srcindices)
-    t_end = time.time()
-    if visualize:
-        print('Time: {:.5f}s'.format(t_end - t_start))
-
-    srcindices = srcindices.get(queue)
-    srcindices_ = srcindices_.get(queue)
-
-    sources = resampler.from_discr.nodes().get(queue)
-    sources_ = resampler.to_discr.nodes().get(queue)
-
-    for i in range(srcindices.nblocks):
-        isrc = srcindices.block_indices(i)
-        isrc_ = srcindices_.block_indices(i)
-
-        for j in range(ndim):
-            assert np.min(sources_[j][isrc_]) <= np.min(sources[j][isrc])
-            assert np.max(sources_[j][isrc_]) >= np.max(sources[j][isrc])
-
-    if visualize:
-        discr = resampler.to_discr
-        _plot_partition_indices(queue, discr, srcindices_,
-                method="elements", use_tree=use_tree, pid="stage2")
-
-
 @pytest.mark.parametrize("ndim", [2, 3])
 @pytest.mark.parametrize("factor", [1.0, 0.6])
 def test_proxy_generator(ctx_factory, ndim, factor, visualize=False):
@@ -253,7 +187,7 @@ def test_proxy_generator(ctx_factory, ndim, factor, visualize=False):
 
     qbx = _build_qbx_discr(queue, ndim=ndim)
     srcindices = _build_block_index(qbx.density_discr,
-            method='nodes', factor=factor)
+            factor=factor)
 
     from pytential.linalg.proxy import ProxyGenerator
     generator = ProxyGenerator(qbx, ratio=1.1)
@@ -274,14 +208,13 @@ def test_proxy_generator(ctx_factory, ndim, factor, visualize=False):
     if visualize:
         if qbx.ambient_dim == 2:
             import matplotlib.pyplot as pt
-            from pytential.qbx.utils import get_centers_on_side
 
             density_nodes = qbx.density_discr.nodes().get(queue)
-            ci = get_centers_on_side(qbx, -1)
+            ci = bind(qbx, sym.expansion_centers(qbx.ambient_dim, -1))(queue)
             ci = np.vstack([c.get(queue) for c in ci])
-            ce = get_centers_on_side(qbx, +1)
+            ce = bind(qbx, sym.expansion_centers(qbx.ambient_dim, +1))(queue)
             ce = np.vstack([c.get(queue) for c in ce])
-            r = qbx._expansion_radii("nsources").get(queue)
+            r = bind(qbx, sym.expansion_radii(qbx.ambient_dim))(queue).get()
 
             for i in range(srcindices.nblocks):
                 isrc = srcindices.block_indices(i)
@@ -333,8 +266,6 @@ def test_proxy_generator(ctx_factory, ndim, factor, visualize=False):
 
                 vis = make_visualizer(queue, discr, 10)
                 filename = "test_proxy_generator_{}d_{:04}.vtu".format(ndim, i)
-                if os.path.isfile(filename):
-                    os.remove(filename)
                 vis.write_vtk_file(filename, [])
 
 
@@ -346,7 +277,7 @@ def test_interaction_points(ctx_factory, ndim, factor, visualize=False):
 
     qbx = _build_qbx_discr(queue, ndim=ndim)
     srcindices = _build_block_index(qbx.density_discr,
-            method='nodes', factor=factor)
+            factor=factor)
 
     # generate proxy points
     from pytential.linalg.proxy import ProxyGenerator
@@ -407,8 +338,6 @@ def test_interaction_points(ctx_factory, ndim, factor, visualize=False):
                 isrc = srcindices.block_indices(i)
                 inbr = nbrindices.block_indices(i)
 
-                # TODO: some way to turn off some of the interpolations
-                # would help visualize this better.
                 marker.fill(0.0)
                 marker[srcindices.indices] = 0.0
                 marker[isrc] = -42.0
@@ -417,9 +346,6 @@ def test_interaction_points(ctx_factory, ndim, factor, visualize=False):
 
                 vis = make_visualizer(queue, qbx.density_discr, 10)
                 filename = "test_area_query_{}d_{:04}.vtu".format(ndim, i)
-                if os.path.isfile(filename):
-                    os.remove(filename)
-
                 vis.write_vtk_file(filename, [
                     ("marker", marker_dev),
                     ])
diff --git a/test/test_matrix.py b/test/test_matrix.py
index d912b2e75f1e932c15e697b130ebfcada5808197..66441657cf0b39a0256d84500fb43120e7362486 100644
--- a/test/test_matrix.py
+++ b/test/test_matrix.py
@@ -40,8 +40,6 @@ from meshmode.mesh.generation import (  # noqa
         ellipse, NArmedStarfish, make_curve_mesh, generate_torus)
 
 from pytential import bind, sym
-from pytential.symbolic.primitives import DEFAULT_SOURCE, DEFAULT_TARGET
-from pytential.symbolic.primitives import QBXSourceStage1, QBXSourceQuadStage2
 
 import pytest
 from pyopencl.tools import (  # noqa
@@ -124,6 +122,8 @@ def _build_block_index(discr, nblks=10, factor=1.0):
 def _build_op(lpot_id,
               k=0,
               ndim=2,
+              source=sym.DEFAULT_SOURCE,
+              target=sym.DEFAULT_TARGET,
               qbx_forced_limit="avg"):
 
     from sumpy.kernel import LaplaceKernel, HelmholtzKernel
@@ -134,7 +134,10 @@ def _build_op(lpot_id,
         knl = LaplaceKernel(ndim)
         knl_kwargs = {}
 
-    lpot_kwargs = {"qbx_forced_limit": qbx_forced_limit}
+    lpot_kwargs = {
+            "qbx_forced_limit": qbx_forced_limit,
+            "source": source,
+            "target": target}
     lpot_kwargs.update(knl_kwargs)
     if lpot_id == 1:
         # scalar single-layer potential
@@ -280,23 +283,33 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id,
     from sympy.core.cache import clear_cache
     clear_cache()
 
+    place_ids = (
+            sym.DOFDescriptor(
+                geometry=sym.DEFAULT_SOURCE,
+                discr_stage=sym.QBX_SOURCE_STAGE1),
+            sym.DOFDescriptor(
+                geometry=sym.DEFAULT_TARGET,
+                discr_stage=sym.QBX_SOURCE_STAGE1),
+            )
     target_order = 2 if ndim == 3 else 7
+
     qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim)
-    op, u_sym, _ = _build_op(lpot_id, ndim=ndim)
+    op, u_sym, _ = _build_op(lpot_id, ndim=ndim,
+            source=place_ids[0],
+            target=place_ids[1])
     index_set = _build_block_index(qbx.density_discr, factor=factor)
 
     from pytential.symbolic.execution import GeometryCollection
-    from pytential.symbolic.execution import _prepare_expr, _prepare_domains
-    places = GeometryCollection(qbx)
+    from pytential.symbolic.execution import _prepare_expr
+    places = GeometryCollection(qbx, auto_where=place_ids)
     expr = _prepare_expr(places, op)
-    domains = _prepare_domains(1, places, None, DEFAULT_SOURCE)
 
     from pytential.symbolic.matrix import P2PMatrixBuilder
     mbuilder = P2PMatrixBuilder(queue,
             dep_expr=u_sym,
             other_dep_exprs=[],
-            dep_source=places[domains[0]],
-            dep_discr=places.get_discretization(domains[0]),
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
             places=places,
             context={},
             exclude_self=True)
@@ -306,8 +319,8 @@ def test_p2p_block_builder(ctx_factory, factor, ndim, lpot_id,
     mbuilder = FarFieldBlockBuilder(queue,
             dep_expr=u_sym,
             other_dep_exprs=[],
-            dep_source=places[domains[0]],
-            dep_discr=places.get_discretization(domains[0]),
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
             places=places,
             index_set=index_set,
             context={},
@@ -349,27 +362,34 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id,
     from sympy.core.cache import clear_cache
     clear_cache()
 
+    place_ids = (
+            sym.DOFDescriptor(
+                geometry=sym.DEFAULT_SOURCE,
+                discr_stage=sym.QBX_SOURCE_STAGE2),
+            sym.DOFDescriptor(
+                geometry=sym.DEFAULT_TARGET,
+                discr_stage=sym.QBX_SOURCE_STAGE2),
+            )
     target_order = 2 if ndim == 3 else 7
-    qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim)
-    op, u_sym, _ = _build_op(lpot_id, ndim=ndim)
 
-    # NOTE: NearFieldBlockBuilder only does stage1/stage1 or stage2/stage2,
-    # so we need to hardcode the discr for MatrixBuilder too, since the
-    # defaults are different
-    where = (QBXSourceStage1(DEFAULT_SOURCE), QBXSourceStage1(DEFAULT_TARGET))
+    qbx = _build_qbx_discr(queue, target_order=target_order, ndim=ndim)
+    op, u_sym, _ = _build_op(lpot_id, ndim=ndim,
+            source=place_ids[0],
+            target=place_ids[1],
+            qbx_forced_limit="avg")
 
     from pytential.symbolic.execution import GeometryCollection, _prepare_expr
-    places = GeometryCollection(qbx, auto_where=where)
+    places = GeometryCollection(qbx, auto_where=place_ids)
     expr = _prepare_expr(places, op)
-    density_discr = places.get_discretization(where[0])
+    density_discr = places.get_discretization(place_ids[0])
     index_set = _build_block_index(density_discr, factor=factor)
 
     from pytential.symbolic.matrix import NearFieldBlockBuilder
     mbuilder = NearFieldBlockBuilder(queue,
             dep_expr=u_sym,
             other_dep_exprs=[],
-            dep_source=places[where[0]],
-            dep_discr=places.get_discretization(where[0]),
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
             places=places,
             index_set=index_set,
             context={})
@@ -379,8 +399,8 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id,
     mbuilder = MatrixBuilder(queue,
             dep_expr=u_sym,
             other_dep_exprs=[],
-            dep_source=places[where[0]],
-            dep_discr=places.get_discretization(where[0]),
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
             places=places,
             context={})
     mat = mbuilder(expr)
@@ -408,13 +428,11 @@ def test_qbx_block_builder(ctx_factory, factor, ndim, lpot_id,
     assert _max_block_error(mat, blk, index_set) < 1.0e-14
 
 
-@pytest.mark.parametrize('place_id',
-        [(DEFAULT_SOURCE, DEFAULT_TARGET),
-         (QBXSourceStage1(DEFAULT_SOURCE),
-          QBXSourceStage1(DEFAULT_TARGET)),
-         (QBXSourceQuadStage2(DEFAULT_SOURCE),
-          QBXSourceQuadStage2(DEFAULT_TARGET))])
-def test_build_matrix_places(ctx_factory, place_id, visualize=False):
+@pytest.mark.parametrize(('source_discr_stage', 'target_discr_stage'),
+        [(sym.QBX_SOURCE_STAGE1, sym.QBX_SOURCE_STAGE1),
+         (sym.QBX_SOURCE_QUAD_STAGE2, sym.QBX_SOURCE_QUAD_STAGE2)])
+def test_build_matrix_places(ctx_factory,
+        source_discr_stage, target_discr_stage, visualize=False):
     ctx = ctx_factory()
     queue = cl.CommandQueue(ctx)
 
@@ -422,38 +440,53 @@ def test_build_matrix_places(ctx_factory, place_id, visualize=False):
     from sympy.core.cache import clear_cache
     clear_cache()
 
+    qbx_forced_limit = -1
+    place_ids = (
+            sym.DOFDescriptor(
+                geometry=sym.DEFAULT_SOURCE,
+                discr_stage=source_discr_stage),
+            sym.DOFDescriptor(
+                geometry=sym.DEFAULT_TARGET,
+                discr_stage=target_discr_stage),
+            )
+
+    # build test operators
     qbx = _build_qbx_discr(queue, nelements=8, target_order=2, ndim=2,
                            curve_f=partial(ellipse, 1.0))
 
-    qbx_forced_limit = -1
     op, u_sym, _ = _build_op(lpot_id=1, ndim=2,
+            source=place_ids[0],
+            target=place_ids[1],
             qbx_forced_limit=qbx_forced_limit)
 
     from pytential.symbolic.execution import GeometryCollection
-    places = GeometryCollection(qbx, auto_where=place_id)
-    source_discr = places.get_discretization(place_id[0])
-    target_discr = places.get_discretization(place_id[1])
+    places = GeometryCollection(qbx, auto_where=place_ids)
+    source_discr = places.get_discretization(place_ids[0])
+    target_discr = places.get_discretization(place_ids[1])
 
     index_set = _build_block_index(source_discr, factor=0.6)
 
-    # build full QBX matrix
-    from pytential.symbolic.execution import build_matrix
-    qbx_mat = build_matrix(queue, qbx, op, u_sym,
-                           auto_where=place_id, domains=place_id[0])
-    qbx_mat = qbx_mat.get(queue)
-
-    assert qbx_mat.shape == (target_discr.nnodes, source_discr.nnodes)
-
-    # build full p2p matrix
     from pytential.symbolic.execution import _prepare_expr
     op = _prepare_expr(places, op)
 
+    # build full QBX matrix
+    from pytential.symbolic.matrix import MatrixBuilder
+    mbuilder = MatrixBuilder(queue,
+            dep_expr=u_sym,
+            other_dep_exprs=[],
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
+            places=places,
+            context={})
+    qbx_mat = mbuilder(op)
+
+    # build full p2p matrix
     from pytential.symbolic.matrix import P2PMatrixBuilder
     mbuilder = P2PMatrixBuilder(queue,
             dep_expr=u_sym,
             other_dep_exprs=[],
-            dep_source=places[place_id[0]],
-            dep_discr=places.get_discretization(place_id[0]),
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
             places=places,
             context={})
     p2p_mat = mbuilder(op)
@@ -465,21 +498,21 @@ def test_build_matrix_places(ctx_factory, place_id, visualize=False):
     mbuilder = NearFieldBlockBuilder(queue,
             dep_expr=u_sym,
             other_dep_exprs=[],
-            dep_source=places[place_id[0]],
-            dep_discr=places.get_discretization(place_id[0]),
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
             places=places,
             index_set=index_set,
             context={})
     mat = mbuilder(op)
-    if place_id[0] is not DEFAULT_SOURCE:
+    if place_ids[0].discr_stage is not None:
         assert _max_block_error(qbx_mat, mat, index_set.get(queue)) < 1.0e-14
 
     from pytential.symbolic.matrix import FarFieldBlockBuilder
     mbuilder = FarFieldBlockBuilder(queue,
             dep_expr=u_sym,
             other_dep_exprs=[],
-            dep_source=places[place_id[0]],
-            dep_discr=places.get_discretization(place_id[0]),
+            dep_source=places.get_geometry(place_ids[0]),
+            dep_discr=places.get_discretization(place_ids[0]),
             places=places,
             index_set=index_set,
             context={},
diff --git a/test/test_maxwell.py b/test/test_maxwell.py
index 69bdc81e19193d3e62a932c20e7e362db96fea93..5adee9d328edd080755bd8e2a206e478b4980d63 100644
--- a/test/test_maxwell.py
+++ b/test/test_maxwell.py
@@ -310,7 +310,7 @@ def test_pec_mfie_extinction(ctx_factory, case, visualize=False):
                     case.fmm_tolerance),
                 fmm_backend=case.fmm_backend
                 ).with_refinement(_expansion_disturbance_tolerance=0.05)
-        h_max = qbx.h_max
+        h_max = bind(qbx, sym.h_max(qbx.ambient_dim))(queue)
 
         scat_discr = qbx.density_discr
         obs_discr = Discretization(
diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py
index 811a912de0ebdf5cc1813ff6e063b1bf1a506d51..2b3dbfa3705bb05813613a19eac3d7d72b4a64af 100644
--- a/test/test_scalar_int_eq.py
+++ b/test/test_scalar_int_eq.py
@@ -622,7 +622,7 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize):
         bc = bind(
                 (point_source, density_discr),
                 sym.normal_derivative(
-                    qbx.ambient_dim, pot_src, where=sym.DEFAULT_TARGET)
+                    qbx.ambient_dim, pot_src, dofdesc=sym.DEFAULT_TARGET)
                 )(queue, charges=source_charges_dev, **concrete_knl_kwargs)
 
     # }}}
@@ -864,8 +864,9 @@ def run_int_eq_test(cl_ctx, queue, case, resolution, visualize):
     class Result(Record):
         pass
 
+    h_max = bind(qbx, sym.h_max(qbx.ambient_dim))(queue)
     return Result(
-            h_max=qbx.h_max,
+            h_max=h_max,
             rel_err_2=rel_err_2,
             rel_err_inf=rel_err_inf,
             rel_td_err_inf=rel_td_err_inf,
diff --git a/test/test_stokes.py b/test/test_stokes.py
index d8d5c821b775ccc3e8057d7f399cfc79452e5bd9..45238e74074a8aa02f25afbd33380bab7cd17fac 100644
--- a/test/test_stokes.py
+++ b/test/test_stokes.py
@@ -268,7 +268,8 @@ def run_exterior_stokes_2d(ctx_factory, nelements,
 
     # }}}
 
-    return qbx.h_max, l2_err
+    h_max = bind(qbx, sym.h_max(qbx.ambient_dim))(queue)
+    return h_max, l2_err
 
 
 @pytest.mark.slowtest
diff --git a/test/test_symbolic.py b/test/test_symbolic.py
index d709f7f0e57cebe40bdf38be90e48e39c8403e5e..40e8382fc8ff5b08013e1ed21726d4f7684f4df5 100644
--- a/test/test_symbolic.py
+++ b/test/test_symbolic.py
@@ -201,13 +201,14 @@ def test_expr_pickling():
 # }}}
 
 
-@pytest.mark.parametrize("source", [
-    sym.DEFAULT_SOURCE,
-    sym.QBXSourceStage1(sym.DEFAULT_SOURCE),
-    sym.QBXSourceStage2(sym.DEFAULT_SOURCE),
-    sym.QBXSourceQuadStage2(sym.DEFAULT_SOURCE)
+@pytest.mark.parametrize(("name", "source_discr_stage", "target_granularity"), [
+    ("default", None, None),
+    ("default-explicit", sym.QBX_SOURCE_STAGE1, sym.GRANULARITY_NODE),
+    ("stage2", sym.QBX_SOURCE_STAGE2, sym.GRANULARITY_NODE),
+    ("stage2-center", sym.QBX_SOURCE_STAGE2, sym.GRANULARITY_CENTER),
+    ("quad", sym.QBX_SOURCE_QUAD_STAGE2, sym.GRANULARITY_NODE)
     ])
-def test_interpolation(ctx_factory, source):
+def test_interpolation(ctx_factory, name, source_discr_stage, target_granularity):
     ctx = ctx_factory()
     queue = cl.CommandQueue(ctx)
 
@@ -227,25 +228,39 @@ def test_interpolation(ctx_factory, source):
             qbx_order=qbx_order,
             fmm_order=False).with_refinement()
 
-    target = sym.QBXSourceQuadStage2(sym.DEFAULT_SOURCE)
+    where = 'test-interpolation'
+    from_dd = sym.DOFDescriptor(
+            geometry=where,
+            discr_stage=source_discr_stage,
+            granularity=sym.GRANULARITY_NODE)
+    to_dd = sym.DOFDescriptor(
+            geometry=where,
+            discr_stage=sym.QBX_SOURCE_QUAD_STAGE2,
+            granularity=target_granularity)
+
     sigma_sym = sym.var("sigma")
-    op_sym = sym.sin(sym.Interpolation(source, target, sigma_sym))
-    bound_op = bind(qbx, op_sym, auto_where=(source, sym.DEFAULT_TARGET))
-
-    quad2_nodes = qbx.quad_stage2_density_discr.nodes().get(queue)
-    if isinstance(source, sym.QBXSourceStage2):
-        nodes = qbx.stage2_density_discr.nodes().get(queue)
-    elif isinstance(source, sym.QBXSourceQuadStage2):
-        nodes = quad2_nodes
+    op_sym = sym.sin(sym.interp(from_dd, to_dd, sigma_sym))
+    bound_op = bind(qbx, op_sym, auto_where=where)
+
+    target_nodes = qbx.quad_stage2_density_discr.nodes().get(queue)
+    if source_discr_stage == sym.QBX_SOURCE_STAGE2:
+        source_nodes = qbx.stage2_density_discr.nodes().get(queue)
+    elif source_discr_stage == sym.QBX_SOURCE_QUAD_STAGE2:
+        source_nodes = target_nodes
     else:
-        nodes = qbx.density_discr.nodes().get(queue)
+        source_nodes = qbx.density_discr.nodes().get(queue)
 
-    sigma_dev = cl.array.to_device(queue, la.norm(nodes, axis=0))
-    sigma_quad2 = np.sin(la.norm(quad2_nodes, axis=0))
-    sigma_quad2_interp = bound_op(queue, sigma=sigma_dev).get(queue)
+    sigma_dev = cl.array.to_device(queue, la.norm(source_nodes, axis=0))
+    sigma_target = np.sin(la.norm(target_nodes, axis=0))
+    sigma_target_interp = bound_op(queue, sigma=sigma_dev).get(queue)
 
-    error = la.norm(sigma_quad2_interp - sigma_quad2) / la.norm(sigma_quad2)
-    assert error < 1.0e-10
+    if name in ('default', 'default-explicit', 'stage2', 'quad'):
+        error = la.norm(sigma_target_interp - sigma_target) / la.norm(sigma_target)
+        assert error < 1.0e-10
+    elif name in ('stage2-center',):
+        assert len(sigma_target_interp) == 2 * len(sigma_target)
+    else:
+        raise ValueError('unknown test case name: {}'.format(name))
 
 
 # You can test individual routines by typing