diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index afb532af7b291ba3e9fffe61de3e9b973e0a2c06..4e03ceca60980f5bf0fccc01df3e9f01df127b5f 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -373,8 +373,17 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
 
     @memoize_method
     def with_refinement(self, target_order=None, kernel_length_scale=None,
-            maxiter=None, visualize=False, _expansion_disturbance_tolerance=None):
+            maxiter=None, visualize=False, refiner=None,
+            _expansion_disturbance_tolerance=None,
+            _force_stage2_uniform_refinement_rounds=None,
+            _scaled_max_curvature_threshold=None):
         """
+        :arg refiner: If the mesh underlying :attr:`density_discr`
+            is itself the result of refinement, then its
+            :class:`meshmode.refinement.Refiner` instance may need to
+            be reused for continued refinement. This argument
+            provides the opportunity to pass in an existing refiner
+            that should be used for continued refinement.
         :returns: a tuple ``(lpot_src, cnx)``, where ``lpot_src`` is a
             :class:`QBXLayerPotentialSource` and ``cnx`` is a
             :class:`meshmode.discretization.connection.DiscretizationConnection`
@@ -395,7 +404,12 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
                     InterpolatoryQuadratureSimplexGroupFactory(target_order),
                     kernel_length_scale=kernel_length_scale,
                     maxiter=maxiter, visualize=visualize,
-                    expansion_disturbance_tolerance=_expansion_disturbance_tolerance)
+                    expansion_disturbance_tolerance=_expansion_disturbance_tolerance,
+                    force_stage2_uniform_refinement_rounds=(
+                        _force_stage2_uniform_refinement_rounds),
+                    scaled_max_curvature_threshold=(
+                        _scaled_max_curvature_threshold),
+                    refiner=refiner)
 
         return lpot, connection
 
@@ -403,8 +417,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
     @memoize_method
     def h_max(self):
         with cl.CommandQueue(self.cl_context) as queue:
-            panel_sizes = self._panel_sizes("npanels").with_queue(queue)
-            return np.asscalar(cl.array.max(panel_sizes).get())
+            quad_res = self._coarsest_quad_resolution("npanels").with_queue(queue)
+            return np.asscalar(cl.array.max(quad_res).get())
 
     # {{{ internal API
 
@@ -414,41 +428,98 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
         return utils.element_centers_of_mass(self.density_discr)
 
     @memoize_method
-    def _fine_panel_centers_of_mass(self):
+    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):
-        if last_dim_length == "npanels":
-            raise ValueError(
-                    "Passing 'npanels' as last_dim_length to _expansion_radii is "
-                    "not allowed. Allowed values are 'nsources' and 'ncenters'.")
-
         with cl.CommandQueue(self.cl_context) as queue:
-                return (self._panel_sizes(last_dim_length).with_queue(queue) * 0.5
-                        ).with_queue(None)
+                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._panel_sizes(last_dim_length).with_queue(queue) * 0.5
+                return (
+                        self._expansion_radii(last_dim_length).with_queue(queue)
+                        * 0.5
                         ).with_queue(None)
 
     @memoize_method
-    def _panel_sizes(self, last_dim_length="npanels"):
+    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
-        return utils.panel_sizes(self.density_discr, last_dim_length)
+        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 _fine_panel_sizes(self, last_dim_length="npanels"):
+    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
-        return utils.panel_sizes(self.stage2_density_discr, last_dim_length)
+        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):
@@ -594,6 +665,8 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
 
         geo_data = self.qbx_fmm_geometry_data(target_discrs_and_qbx_sides)
 
+        # geo_data.plot()
+
         # FIXME Exert more positive control over geo_data attribute lifetimes using
         # geo_data.<method>.clear_cache(geo_data).
 
diff --git a/pytential/qbx/fmm.py b/pytential/qbx/fmm.py
index 037f818893a77edbad9a032a3868ac0dec8af514..d3622e5923e2846acd02285303e8e76d7054ad02 100644
--- a/pytential/qbx/fmm.py
+++ b/pytential/qbx/fmm.py
@@ -33,6 +33,8 @@ from sumpy.fmm import (SumpyExpansionWranglerCodeContainer,
 from pytools import memoize_method
 from pytential.qbx.interactions import P2QBXLFromCSR, M2QBXL, L2QBXL, QBXL2P
 
+from pytools import log_process, ProcessLogger
+
 import logging
 logger = logging.getLogger(__name__)
 
@@ -192,6 +194,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
     # {{{ qbx-related
 
+    @log_process(logger)
     def form_global_qbx_locals(self, src_weights):
         local_exps = self.qbx_local_expansion_zeros()
 
@@ -228,6 +231,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
         return result
 
+    @log_process(logger)
     def translate_box_multipoles_to_qbx_local(self, multipole_exps):
         qbx_expansions = self.qbx_local_expansion_zeros()
 
@@ -276,6 +280,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
         return qbx_expansions
 
+    @log_process(logger)
     def translate_box_local_to_qbx_local(self, local_exps):
         qbx_expansions = self.qbx_local_expansion_zeros()
 
@@ -320,6 +325,7 @@ QBXFMMGeometryData.non_qbx_box_target_lists`),
 
         return qbx_expansions
 
+    @log_process(logger)
     def eval_qbx_expansions(self, qbx_expansions):
         pot = self.full_output_zeros()
 
@@ -381,16 +387,12 @@ def drive_fmm(expansion_wrangler, src_weights):
     # Interface guidelines: Attributes of the tree are assumed to be known
     # to the expansion wrangler and should not be passed.
 
-    from time import time
-    start_time = time()
-    logger.info("start qbx fmm")
+    fmm_proc = ProcessLogger(logger, "qbx fmm")
 
-    logger.info("reorder source weights")
     src_weights = wrangler.reorder_sources(src_weights)
 
     # {{{ construct local multipoles
 
-    logger.info("construct local multipoles")
     mpole_exps = wrangler.form_multipoles(
             traversal.level_start_source_box_nrs,
             traversal.source_boxes,
@@ -400,7 +402,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ propagate multipoles upward
 
-    logger.info("propagate multipoles upward")
     wrangler.coarsen_multipoles(
             traversal.level_start_source_parent_box_nrs,
             traversal.source_parent_boxes,
@@ -410,7 +411,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ direct evaluation from neighbor source boxes ("list 1")
 
-    logger.info("direct evaluation from neighbor source boxes ('list 1')")
     non_qbx_potentials = wrangler.eval_direct(
             traversal.target_boxes,
             traversal.neighbor_source_boxes_starts,
@@ -421,7 +421,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ translate separated siblings' ("list 2") mpoles to local
 
-    logger.info("translate separated siblings' ('list 2') mpoles to local")
     local_exps = wrangler.multipole_to_local(
             traversal.level_start_target_or_target_parent_box_nrs,
             traversal.target_or_target_parent_boxes,
@@ -433,8 +432,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ evaluate sep. smaller mpoles ("list 3") at particles
 
-    logger.info("evaluate sep. smaller mpoles at particles ('list 3 far')")
-
     # (the point of aiming this stage at particles is specifically to keep its
     # contribution *out* of the downward-propagating local expansions)
 
@@ -450,8 +447,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ form locals for separated bigger source boxes ("list 4")
 
-    logger.info("form locals for separated bigger source boxes ('list 4 far')")
-
     local_exps = local_exps + wrangler.form_locals(
             traversal.level_start_target_or_target_parent_box_nrs,
             traversal.target_or_target_parent_boxes,
@@ -466,7 +461,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ propagate local_exps downward
 
-    logger.info("propagate local_exps downward")
     wrangler.refine_locals(
             traversal.level_start_target_or_target_parent_box_nrs,
             traversal.target_or_target_parent_boxes,
@@ -476,7 +470,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ evaluate locals
 
-    logger.info("evaluate locals")
     non_qbx_potentials = non_qbx_potentials + wrangler.eval_locals(
             traversal.level_start_target_box_nrs,
             traversal.target_boxes,
@@ -486,19 +479,14 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ wrangle qbx expansions
 
-    logger.info("form global qbx expansions from list 1")
     qbx_expansions = wrangler.form_global_qbx_locals(src_weights)
 
-    logger.info("translate from list 3 multipoles to qbx local expansions")
     qbx_expansions = qbx_expansions + \
             wrangler.translate_box_multipoles_to_qbx_local(mpole_exps)
 
-    logger.info("translate from box local expansions to contained "
-            "qbx local expansions")
     qbx_expansions = qbx_expansions + \
             wrangler.translate_box_local_to_qbx_local(local_exps)
 
-    logger.info("evaluate qbx local expansions")
     qbx_potentials = wrangler.eval_qbx_expansions(
             qbx_expansions)
 
@@ -506,8 +494,6 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # {{{ reorder potentials
 
-    logger.info("reorder potentials")
-
     nqbtl = geo_data.non_qbx_box_target_lists()
 
     all_potentials_in_tree_order = wrangler.full_output_zeros()
@@ -528,7 +514,7 @@ def drive_fmm(expansion_wrangler, src_weights):
 
     # }}}
 
-    logger.info("qbx fmm complete in %.2f s" % (time() - start_time))
+    fmm_proc.done()
 
     return result
 
diff --git a/pytential/qbx/fmmlib.py b/pytential/qbx/fmmlib.py
index 578dadce208da52b81f1e5014381e0aa04df4a17..2b6cbf88ce53364e4b7efc0b1ad01529eaaea0c9 100644
--- a/pytential/qbx/fmmlib.py
+++ b/pytential/qbx/fmmlib.py
@@ -30,6 +30,8 @@ from boxtree.pyfmmlib_integration import FMMLibExpansionWrangler
 from sumpy.kernel import LaplaceKernel, HelmholtzKernel
 
 
+from pytools import log_process
+
 import logging
 logger = logging.getLogger(__name__)
 
@@ -298,6 +300,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler):
 
     # {{{ p2qbxl
 
+    @log_process(logger)
     def form_global_qbx_locals(self, src_weights):
         geo_data = self.geo_data
         trav = geo_data.traversal()
@@ -348,6 +351,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler):
 
     # {{{ m2qbxl
 
+    @log_process(logger)
     def translate_box_multipoles_to_qbx_local(self, multipole_exps):
         qbx_exps = self.qbx_local_expansion_zeros()
 
@@ -458,6 +462,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler):
 
     # }}}
 
+    @log_process(logger)
     def translate_box_local_to_qbx_local(self, local_exps):
         qbx_expansions = self.qbx_local_expansion_zeros()
 
@@ -518,6 +523,7 @@ class QBXFMMLibExpansionWrangler(FMMLibExpansionWrangler):
 
         return qbx_expansions
 
+    @log_process(logger)
     def eval_qbx_expansions(self, qbx_expansions):
         output = self.full_output_zeros()
 
diff --git a/pytential/qbx/geometry.py b/pytential/qbx/geometry.py
index 47c578b764f39287e1dce2c39ec32d06e3ca38b1..ee4421848fb7ba0388e76682be256727d47313ac 100644
--- a/pytential/qbx/geometry.py
+++ b/pytential/qbx/geometry.py
@@ -32,11 +32,11 @@ from boxtree.tools import DeviceDataRecord
 import loopy as lp
 from loopy.version import MOST_RECENT_LANGUAGE_VERSION
 from cgen import Enum
-from time import time
 
 
 from pytential.qbx.utils import TreeCodeContainerMixin
 
+from pytools import log_process
 
 import logging
 logger = logging.getLogger(__name__)
@@ -670,6 +670,7 @@ class QBXFMMGeometryData(object):
         return result.with_queue(None)
 
     @memoize_method
+    @log_process(logger)
     def global_qbx_centers(self):
         """Build a list of indices of QBX centers that use global QBX.  This
         indexes into the global list of targets, (see :meth:`target_info`) of
@@ -682,9 +683,6 @@ class QBXFMMGeometryData(object):
         user_target_to_center = self.user_target_to_center()
 
         with cl.CommandQueue(self.cl_context) as queue:
-            logger.debug("find global qbx centers: start")
-            center_find_start_time = time()
-
             tgt_assoc_result = (
                     user_target_to_center.with_queue(queue)[self.ncenters:])
 
@@ -709,15 +707,6 @@ class QBXFMMGeometryData(object):
                         ],
                     queue=queue)
 
-            center_find_elapsed = time() - center_find_start_time
-            if center_find_elapsed > 0.1:
-                done_logger = logger.info
-            else:
-                done_logger = logger.debug
-
-            done_logger("find global qbx centers: done after %g seconds",
-                    center_find_elapsed)
-
             if self.debug:
                 logger.debug(
                         "find global qbx centers: using %d/%d centers"
@@ -767,6 +756,7 @@ class QBXFMMGeometryData(object):
         return result.with_queue(None)
 
     @memoize_method
+    @log_process(logger)
     def center_to_tree_targets(self):
         """Return a :class:`CenterToTargetList`. See :meth:`target_to_center`
         for the reverse look-up table with targets in user order.
@@ -777,9 +767,6 @@ class QBXFMMGeometryData(object):
         user_ttc = self.user_target_to_center()
 
         with cl.CommandQueue(self.cl_context) as queue:
-            logger.debug("build center -> targets lookup table: start")
-            c2t_start_time = time()
-
             tree_ttc = cl.array.empty_like(user_ttc).with_queue(queue)
             tree_ttc[self.tree().sorted_target_ids] = user_ttc
 
@@ -802,14 +789,6 @@ class QBXFMMGeometryData(object):
                             filtered_tree_ttc, filtered_target_ids,
                             self.ncenters, tree_ttc.dtype)
 
-            c2t_elapsed = time() - c2t_start_time
-            if c2t_elapsed > 0.1:
-                done_logger = logger.info
-            else:
-                done_logger = logger.debug
-            done_logger("build center -> targets lookup table: "
-                    "done after %g seconds", c2t_elapsed)
-
             result = CenterToTargetList(
                     starts=center_target_starts,
                     lists=targets_sorted_by_center).with_queue(None)
@@ -817,6 +796,7 @@ class QBXFMMGeometryData(object):
             return result
 
     @memoize_method
+    @log_process(logger)
     def non_qbx_box_target_lists(self):
         """Build a list of targets per box that don't need to bother with QBX.
         Returns a :class:`boxtree.tree.FilteredTargetListsInTreeOrder`.
@@ -827,9 +807,6 @@ class QBXFMMGeometryData(object):
         """
 
         with cl.CommandQueue(self.cl_context) as queue:
-            logger.debug("find non-qbx box target lists: start")
-            nonqbx_start_time = time()
-
             flags = (self.user_target_to_center().with_queue(queue)
                     == target_state.NO_QBX_NEEDED)
 
@@ -845,14 +822,6 @@ class QBXFMMGeometryData(object):
             plfilt = self.code_getter.particle_list_filter()
             result = plfilt.filter_target_lists_in_tree_order(queue, tree, flags)
 
-            nonqbx_elapsed = time() - nonqbx_start_time
-            if nonqbx_elapsed > 0.1:
-                done_logger = logger.info
-            else:
-                done_logger = logger.debug
-            done_logger("find non-qbx box target lists: done after %g seconds",
-                    nonqbx_elapsed)
-
             return result.with_queue(None)
 
     # {{{ plotting (for debugging)
@@ -870,12 +839,14 @@ class QBXFMMGeometryData(object):
             This only works for two-dimensional geometries.
         """
 
+        import matplotlib.pyplot as pt
+        pt.clf()
+
         dims = self.tree().targets.shape[0]
         if dims != 2:
             raise ValueError("only 2-dimensional geometry info can be plotted")
 
         with cl.CommandQueue(self.cl_context) as queue:
-            import matplotlib.pyplot as pt
             from meshmode.discretization.visualization import draw_curve
             draw_curve(self.lpot_source.quad_stage2_density_discr)
 
@@ -960,8 +931,10 @@ class QBXFMMGeometryData(object):
             # }}}
 
             pt.gca().set_aspect("equal")
-            pt.legend()
-            pt.show()
+            #pt.legend()
+            pt.savefig(
+                    "geodata-stage2-nelem%d.pdf"
+                    % self.lpot_source.stage2_density_discr.mesh.nelements)
 
     # }}}
 
diff --git a/pytential/qbx/refinement.py b/pytential/qbx/refinement.py
index 08539c6899e696dd8af6778bdefdd41de33d87a8..2ab458a2b12b6d749791b8b23d771a12a627b069 100644
--- a/pytential/qbx/refinement.py
+++ b/pytential/qbx/refinement.py
@@ -39,6 +39,8 @@ from pytential.qbx.utils import (
         QBX_TREE_C_PREAMBLE, QBX_TREE_MAKO_DEFS, TreeWranglerBase,
         TreeCodeContainerMixin)
 
+from pytools import ProcessLogger, log_process
+
 import logging
 logger = logging.getLogger(__name__)
 
@@ -242,23 +244,24 @@ class RefinerCodeContainer(TreeCodeContainerMixin):
                 extra_type_aliases=(("particle_id_t", particle_id_dtype),))
 
     @memoize_method
-    def kernel_length_scale_to_panel_size_ratio_checker(self):
+    def element_prop_threshold_checker(self):
         knl = lp.make_kernel(
-            "{[panel]: 0<=panel<npanels}",
+            "{[ielement]: 0<=ielement<nelements}",
             """
-            for panel
-                <> oversize = panel_sizes[panel] > kernel_length_scale
-                if oversize
-                    refine_flags[panel] = 1
+            for ielement
+                <> over_threshold = element_property[ielement] > threshold
+                if over_threshold
+                    refine_flags[ielement] = 1
                     refine_flags_updated = 1 {id=write_refine_flags_updated}
                 end
             end
             """,
             options="return_dict",
             silenced_warnings="write_race(write_refine_flags_updated)",
-            name="refine_kernel_length_scale_to_panel_size_ratio",
+            name="refine_kernel_length_scale_to_quad_resolution_ratio",
             lang_version=MOST_RECENT_LANGUAGE_VERSION)
-        knl = lp.split_iname(knl, "panel", 128, inner_tag="l.0", outer_tag="g.0")
+
+        knl = lp.split_iname(knl, "ielement", 128, inner_tag="l.0", outer_tag="g.0")
         return knl
 
     def get_wrangler(self, queue):
@@ -280,6 +283,7 @@ class RefinerWrangler(TreeWranglerBase):
 
     # {{{ check subroutines for conditions 1-3
 
+    @log_process(logger)
     def check_expansion_disks_undisturbed_by_sources(self,
             lpot_source, tree, peer_lists,
             expansion_disturbance_tolerance,
@@ -298,9 +302,6 @@ class RefinerWrangler(TreeWranglerBase):
                 tree.particle_id_dtype,
                 max_levels)
 
-        logger.info("refiner: checking that expansion disk is "
-                "undisturbed by sources")
-
         if debug:
             npanels_to_refine_prev = cl.array.sum(refine_flags).get()
 
@@ -338,10 +339,9 @@ class RefinerWrangler(TreeWranglerBase):
                 logger.debug("refiner: found {} panel(s) to refine".format(
                     npanels_to_refine - npanels_to_refine_prev))
 
-        logger.info("refiner: done checking center is closest to orig panel")
-
         return found_panel_to_refine.get()[0] == 1
 
+    @log_process(logger)
     def check_sufficient_source_quadrature_resolution(
             self, lpot_source, tree, peer_lists, refine_flags, debug,
             wait_for=None):
@@ -360,14 +360,11 @@ class RefinerWrangler(TreeWranglerBase):
         if debug:
             npanels_to_refine_prev = cl.array.sum(refine_flags).get()
 
-        logger.info("refiner: checking sufficient quadrature resolution")
-
         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._fine_panel_sizes("npanels")
-                .with_queue(self.queue) / 4)
+        source_danger_zone_radii_by_panel = \
+                lpot_source._source_danger_zone_radii("npanels")
         unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
 
         evt = knl(
@@ -396,24 +393,21 @@ class RefinerWrangler(TreeWranglerBase):
                 logger.debug("refiner: found {} panel(s) to refine".format(
                     npanels_to_refine - npanels_to_refine_prev))
 
-        logger.info("refiner: done checking sufficient quadrature resolution")
-
         return found_panel_to_refine.get()[0] == 1
 
-    def check_kernel_length_scale_to_panel_size_ratio(self, lpot_source,
-            kernel_length_scale, refine_flags, debug, wait_for=None):
-        knl = self.code_container.kernel_length_scale_to_panel_size_ratio_checker()
-
-        logger.info("refiner: checking kernel length scale to panel size ratio")
+    def check_element_prop_threshold(self, element_property, threshold, refine_flags,
+            debug, wait_for=None):
+        knl = self.code_container.element_prop_threshold_checker()
 
         if debug:
             npanels_to_refine_prev = cl.array.sum(refine_flags).get()
 
         evt, out = knl(self.queue,
-                       panel_sizes=lpot_source._panel_sizes("npanels"),
+                       element_property=element_property,
+                       # lpot_source._coarsest_quad_resolution("npanels")),
                        refine_flags=refine_flags,
                        refine_flags_updated=np.array(0),
-                       kernel_length_scale=np.array(kernel_length_scale),
+                       threshold=np.array(threshold),
                        wait_for=wait_for)
 
         cl.wait_for_events([evt])
@@ -424,8 +418,6 @@ class RefinerWrangler(TreeWranglerBase):
                 logger.debug("refiner: found {} panel(s) to refine".format(
                     npanels_to_refine - npanels_to_refine_prev))
 
-        logger.info("refiner: done checking kernel length scale to panel size ratio")
-
         return (out["refine_flags_updated"].get() == 1).all()
 
     # }}}
@@ -438,13 +430,10 @@ class RefinerWrangler(TreeWranglerBase):
             refine_flags = refine_flags.get(self.queue)
         refine_flags = refine_flags.astype(np.bool)
 
-        logger.info("refiner: calling meshmode")
-
-        refiner.refine(refine_flags)
-        from meshmode.discretization.connection import make_refinement_connection
-        conn = make_refinement_connection(refiner, density_discr, factory)
-
-        logger.info("refiner: done calling meshmode")
+        with ProcessLogger(logger, "refine mesh"):
+            refiner.refine(refine_flags)
+            from meshmode.discretization.connection import make_refinement_connection
+            conn = make_refinement_connection(refiner, density_discr, factory)
 
         return conn
 
@@ -476,8 +465,11 @@ def make_empty_refine_flags(queue, lpot_source, use_stage2_discr=False):
 
 def refine_for_global_qbx(lpot_source, wrangler,
         group_factory, kernel_length_scale=None,
-        refine_flags=None, debug=None, maxiter=None,
-        visualize=None, expansion_disturbance_tolerance=None):
+        force_stage2_uniform_refinement_rounds=None,
+        scaled_max_curvature_threshold=None,
+        debug=None, maxiter=None,
+        visualize=None, expansion_disturbance_tolerance=None,
+        refiner=None):
     """
     Entry point for calling the refiner.
 
@@ -492,11 +484,6 @@ def refine_for_global_qbx(lpot_source, wrangler,
     :arg kernel_length_scale: The kernel length scale, or *None* if not
         applicable. All panels are refined to below this size.
 
-    :arg refine_flags: A :class:`pyopencl.array.Array` indicating which
-        panels should get refined initially, or `None` if no initial
-        refinement should be done. Should have size equal to the number of
-        panels. See also :func:`make_empty_refine_flags()`.
-
     :arg maxiter: The maximum number of refiner iterations.
 
     :returns: A tuple ``(lpot_source, *conn*)`` where ``lpot_source`` is the
@@ -515,23 +502,24 @@ def refine_for_global_qbx(lpot_source, wrangler,
     if expansion_disturbance_tolerance is None:
         expansion_disturbance_tolerance = 0.025
 
+    if force_stage2_uniform_refinement_rounds is None:
+        force_stage2_uniform_refinement_rounds = 0
+
     # TODO: Stop doing redundant checks by avoiding panels which no longer need
     # refinement.
 
-    from meshmode.mesh.refinement import Refiner
+    from meshmode.mesh.refinement import RefinerWithoutAdjacency
     from meshmode.discretization.connection import (
             ChainedDiscretizationConnection, make_same_mesh_connection)
 
-    refiner = Refiner(lpot_source.density_discr.mesh)
-    connections = []
+    if refiner is not None:
+        assert refiner.get_current_mesh() == lpot_source.density_discr.mesh
+    else:
+        # We may be handed a mesh that's already non-conforming, we don't rely
+        # on adjacency, and the no-adjacency refiner is faster.
+        refiner = RefinerWithoutAdjacency(lpot_source.density_discr.mesh)
 
-    # Do initial refinement.
-    if refine_flags is not None:
-        conn = wrangler.refine(
-                lpot_source.density_discr, refiner, refine_flags, group_factory,
-                debug)
-        connections.append(conn)
-        lpot_source = lpot_source.copy(density_discr=conn.to_discr)
+    connections = []
 
     # {{{ first stage refinement
 
@@ -541,7 +529,7 @@ def refine_for_global_qbx(lpot_source, wrangler,
 
         discr = lpot_source.density_discr
         from meshmode.discretization.visualization import make_visualizer
-        vis = make_visualizer(wrangler.queue, discr, 10)
+        vis = make_visualizer(wrangler.queue, discr, 3)
 
         flags = flags.get().astype(np.bool)
         nodes_flags = np.zeros(discr.nnodes)
@@ -598,32 +586,67 @@ def refine_for_global_qbx(lpot_source, wrangler,
             warn_max_iterations()
             break
 
-        # Build tree and auxiliary data.
-        # FIXME: The tree should not have to be rebuilt at each iteration.
-        tree = wrangler.build_tree(lpot_source)
-        peer_lists = wrangler.find_peer_lists(tree)
         refine_flags = make_empty_refine_flags(wrangler.queue, lpot_source)
 
-        # Check condition 1.
-        has_disturbed_expansions = \
-                wrangler.check_expansion_disks_undisturbed_by_sources(
-                        lpot_source, tree, peer_lists,
-                        expansion_disturbance_tolerance,
-                        refine_flags, debug)
-        if has_disturbed_expansions:
-            iter_violated_criteria.append("disturbed expansions")
-            visualize_refinement(niter, "disturbed-expansions", refine_flags)
-
-        # Check condition 3.
         if kernel_length_scale is not None:
-
-            violates_kernel_length_scale = \
-                    wrangler.check_kernel_length_scale_to_panel_size_ratio(
-                            lpot_source, kernel_length_scale, refine_flags, debug)
-
-            if violates_kernel_length_scale:
-                iter_violated_criteria.append("kernel length scale")
-                visualize_refinement(niter, "kernel-length-scale", refine_flags)
+            with ProcessLogger(logger,
+                    "checking kernel length scale to panel size ratio"):
+
+                violates_kernel_length_scale = \
+                        wrangler.check_element_prop_threshold(
+                                element_property=(
+                                    lpot_source._coarsest_quad_resolution(
+                                        "npanels")),
+                                threshold=kernel_length_scale,
+                                refine_flags=refine_flags, debug=debug)
+
+                if violates_kernel_length_scale:
+                    iter_violated_criteria.append("kernel length scale")
+                    visualize_refinement(niter, "kernel-length-scale", refine_flags)
+
+        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")
+
+                violates_scaled_max_curv = \
+                        wrangler.check_element_prop_threshold(
+                                element_property=scaled_max_curv,
+                                threshold=scaled_max_curvature_threshold,
+                                refine_flags=refine_flags, debug=debug)
+
+                if violates_scaled_max_curv:
+                    iter_violated_criteria.append("curvature")
+                    visualize_refinement(niter, "curvature", refine_flags)
+
+        if not iter_violated_criteria:
+            # Only start building trees once the simple length-based criteria
+            # are happy.
+
+            # Build tree and auxiliary data.
+            # FIXME: The tree should not have to be rebuilt at each iteration.
+            tree = wrangler.build_tree(lpot_source)
+            peer_lists = wrangler.find_peer_lists(tree)
+
+            has_disturbed_expansions = \
+                    wrangler.check_expansion_disks_undisturbed_by_sources(
+                            lpot_source, tree, peer_lists,
+                            expansion_disturbance_tolerance,
+                            refine_flags, debug)
+            if has_disturbed_expansions:
+                iter_violated_criteria.append("disturbed expansions")
+                visualize_refinement(niter, "disturbed-expansions", refine_flags)
+
+            del tree
+            del peer_lists
 
         if iter_violated_criteria:
             violated_criteria.append(" and ".join(iter_violated_criteria))
@@ -634,9 +657,7 @@ def refine_for_global_qbx(lpot_source, wrangler,
             connections.append(conn)
             lpot_source = lpot_source.copy(density_discr=conn.to_discr)
 
-        del tree
         del refine_flags
-        del peer_lists
 
     # }}}
 
@@ -686,6 +707,18 @@ def refine_for_global_qbx(lpot_source, wrangler,
         del refine_flags
         del peer_lists
 
+    for round in range(force_stage2_uniform_refinement_rounds):
+        conn = wrangler.refine(
+                stage2_density_discr,
+                refiner,
+                np.ones(stage2_density_discr.mesh.nelements, dtype=np.bool),
+                group_factory, debug)
+        stage2_density_discr = conn.to_discr
+        fine_connections.append(conn)
+        lpot_source = lpot_source.copy(
+                to_refined_connection=ChainedDiscretizationConnection(
+                    fine_connections))
+
     # }}}
 
     lpot_source = lpot_source.copy(debug=debug, _refined_for_global_qbx=True)
diff --git a/pytential/qbx/target_assoc.py b/pytential/qbx/target_assoc.py
index 8d85a5dcf6715182ada6918332303dc48d000b9c..f51e920fe1663db4f6d9dee4f0b6429b150eceb1 100644
--- a/pytential/qbx/target_assoc.py
+++ b/pytential/qbx/target_assoc.py
@@ -39,11 +39,11 @@ from cgen import Enum
 from pytential.qbx.utils import (
     QBX_TREE_C_PREAMBLE, QBX_TREE_MAKO_DEFS, TreeWranglerBase,
     TreeCodeContainerMixin)
-from time import time
 
 
 unwrap_args = AreaQueryElementwiseTemplate.unwrap_args
 
+from pytools import log_process
 
 import logging
 logger = logging.getLogger(__name__)
@@ -439,6 +439,7 @@ class TargetAssociationWrangler(TreeWranglerBase):
         self.code_container = code_container
         self.queue = queue
 
+    @log_process(logger)
     def mark_targets(self, tree, peer_lists, lpot_source, target_status,
                      debug, wait_for=None):
         # Round up level count--this gets included in the kernel as
@@ -498,9 +499,6 @@ class TargetAssociationWrangler(TreeWranglerBase):
                 wait_for=wait_for)
         wait_for = [evt]
 
-        logger.debug("target association: marking targets close to panels")
-        mark_start_time = time()
-
         tunnel_radius_by_source = lpot_source._close_target_tunnel_radius("nsources")
 
         evt = knl(
@@ -529,16 +527,9 @@ class TargetAssociationWrangler(TreeWranglerBase):
 
         cl.wait_for_events([evt])
 
-        mark_elapsed = time() - mark_start_time
-        if mark_elapsed > 0.1:
-            done_logger = logger.info
-        else:
-            done_logger = logger.debug
-        done_logger("target association: done marking targets close to panels "
-                "after %g seconds", mark_elapsed)
-
         return (found_target_close_to_panel == 1).all().get()
 
+    @log_process(logger)
     def try_find_centers(self, tree, peer_lists, lpot_source,
                          target_status, target_flags, target_assoc,
                          target_association_tolerance, debug, wait_for=None):
@@ -590,9 +581,6 @@ class TargetAssociationWrangler(TreeWranglerBase):
 
         wait_for.extend(min_dist_to_center.events)
 
-        logger.debug("target association: finding centers for targets")
-        center_find_start_time = time()
-
         evt = knl(
             *unwrap_args(
                 tree, peer_lists,
@@ -623,15 +611,7 @@ class TargetAssociationWrangler(TreeWranglerBase):
 
         cl.wait_for_events([evt])
 
-        center_find_elapsed = time() - center_find_start_time
-        if center_find_elapsed > 0.1:
-            done_logger = logger.info
-        else:
-            done_logger = logger.debug
-
-        done_logger("target association: done finding centers for targets "
-                "after %g seconds", center_find_elapsed)
-
+    @log_process(logger)
     def mark_panels_for_refinement(self, tree, peer_lists, lpot_source,
                                    target_status, refine_flags, debug,
                                    wait_for=None):
@@ -669,9 +649,6 @@ class TargetAssociationWrangler(TreeWranglerBase):
                 wait_for=wait_for)
         wait_for = [evt]
 
-        logger.debug("target association: marking panels for refinement")
-        mark_start_time = time()
-
         evt = knl(
             *unwrap_args(
                 tree, peer_lists,
@@ -701,15 +678,6 @@ class TargetAssociationWrangler(TreeWranglerBase):
 
         cl.wait_for_events([evt])
 
-        mark_elapsed = time() - mark_start_time
-        if mark_elapsed > 0.1:
-            done_logger = logger.info
-        else:
-            done_logger = logger.debug
-
-        done_logger("target association: done marking panels for refinement "
-                "after %g seconds", mark_elapsed)
-
         return (found_panel_to_refine == 1).all().get()
 
     def make_target_flags(self, target_discrs_and_qbx_sides):
diff --git a/pytential/qbx/utils.py b/pytential/qbx/utils.py
index 6673b450bb7ed3cb2cdf043e95963c0930fab89a..b0ffc066035c0c8f1aea690ecc5dcb2618367275 100644
--- a/pytential/qbx/utils.py
+++ b/pytential/qbx/utils.py
@@ -37,6 +37,8 @@ from loopy.version import MOST_RECENT_LANGUAGE_VERSION
 import logging
 logger = logging.getLogger(__name__)
 
+from pytools import log_process
+
 
 # {{{ c and mako snippets
 
@@ -205,68 +207,51 @@ class TreeWranglerBase(object):
 # }}}
 
 
-# {{{ panel sizes
+# {{{ to_last_dim_length
 
-def panel_sizes(discr, last_dim_length):
-    if last_dim_length not in ("nsources", "ncenters", "npanels"):
-        raise ValueError(
-                "invalid value of last_dim_length: %s" % last_dim_length)
+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*.
+    """
 
-    # To get the panel size this does the equivalent of (∫ 1 ds)**(1/dim).
-    # FIXME: Kernel optimizations
+    queue = queue or vec.queue
 
-    if last_dim_length == "nsources" or last_dim_length == "ncenters":
-        knl = lp.make_kernel(
-            "{[i,j,k]: 0<=i<nelements and 0<=j,k<nunit_nodes}",
-            "panel_sizes[i,j] = sum(k, ds[i,k])**(1/dim)",
-            name="compute_size",
-            lang_version=MOST_RECENT_LANGUAGE_VERSION)
+    if last_dim_length == "nsources":
+        return vec
 
-        def panel_size_view(discr, group_nr):
-            return discr.groups[group_nr].view
+    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,j]: 0<=i<nelements and 0<=j<nunit_nodes}",
-            "panel_sizes[i] = sum(j, ds[i,j])**(1/dim)",
-            name="compute_size",
-            lang_version=MOST_RECENT_LANGUAGE_VERSION)
-        from functools import partial
-
-        def panel_size_view(discr, group_nr):
-            return partial(el_view, discr, group_nr)
+            "{[i,k]: 0<=i<nelements}",
+            "result[i] = a[i,0]",
+            [
+                lp.GlobalArg("a", shape="nelements, nunit_nodes", dtype=lp.auto),
+                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")
 
-    knl = lp.fix_parameters(knl, dim=discr.dim)
-
-    with cl.CommandQueue(discr.cl_context) as queue:
-        from pytential import bind, sym
-        ds = bind(
-                discr,
-                sym.area_element(ambient_dim=discr.ambient_dim, dim=discr.dim)
-                * sym.QWeight()
-                )(queue)
-        panel_sizes = cl.array.empty(
-            queue, discr.nnodes
-            if last_dim_length in ("nsources", "ncenters")
-            else discr.mesh.nelements, discr.real_dtype)
-        for group_nr, group in enumerate(discr.groups):
-            _, (result,) = knl(queue,
-                nelements=group.nelements,
-                nunit_nodes=group.nunit_nodes,
-                ds=group.view(ds),
-                panel_sizes=panel_size_view(
-                    discr, group_nr)(panel_sizes))
-        panel_sizes.finish()
-        if last_dim_length == "ncenters":
-            from pytential.qbx.utils import get_interleaver_kernel
-            knl = get_interleaver_kernel(discr.real_dtype)
-            _, (panel_sizes,) = knl(queue, dstlen=2*discr.nnodes,
-                                    src1=panel_sizes, src2=panel_sizes)
-        return panel_sizes.with_queue(None)
-
 # }}}
 
 
@@ -331,11 +316,10 @@ def el_view(discr, group_nr, global_array):
     where *nelements* is the global (per-discretization) node count.
     """
 
-    group = discr.groups[group_nr]
-    el_nr_base = sum(group.nelements for group in discr.groups[:group_nr])
+    group = discr.mesh.groups[group_nr]
 
     return global_array[
-        ..., el_nr_base:el_nr_base + group.nelements] \
+        ..., group.element_nr_base:group.element_nr_base + group.nelements] \
         .reshape(
             global_array.shape[:-1]
             + (group.nelements,))
@@ -461,6 +445,7 @@ class TreeWithQBXMetadata(Tree):
 MAX_REFINE_WEIGHT = 64
 
 
+@log_process(logger)
 def build_tree_with_qbx_metadata(
         queue, tree_builder, particle_list_filter, lpot_source, targets_list=(),
         use_stage2_discr=False):
@@ -490,8 +475,6 @@ def build_tree_with_qbx_metadata(
     # - then panels (=centers of mass)
     # - then targets
 
-    logger.info("start building tree with qbx metadata")
-
     sources = (
             lpot_source.density_discr.nodes()
             if not use_stage2_discr
@@ -502,7 +485,7 @@ def build_tree_with_qbx_metadata(
     centers_of_mass = (
             lpot_source._panel_centers_of_mass()
             if not use_stage2_discr
-            else lpot_source._fine_panel_centers_of_mass())
+            else lpot_source._stage2_panel_centers_of_mass())
 
     targets = (tgt.nodes() for tgt in targets_list)
 
@@ -603,8 +586,6 @@ def build_tree_with_qbx_metadata(
 
     tree_attrs.update(particle_classes)
 
-    logger.info("done building tree with qbx metadata")
-
     return TreeWithQBXMetadata(
         qbx_panel_to_source_starts=qbx_panel_to_source_starts,
         qbx_panel_to_center_starts=qbx_panel_to_center_starts,
diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py
index 96cac5c4f5074bfb74f8a7460855894210095c26..c0fe01535c047758d9240b0af306e44e5ee12506 100644
--- a/pytential/symbolic/execution.py
+++ b/pytential/symbolic/execution.py
@@ -33,6 +33,10 @@ import pyopencl as cl
 import pyopencl.array  # noqa
 import pyopencl.clmath  # noqa
 
+from loopy.version import MOST_RECENT_LANGUAGE_VERSION
+
+from pytools import memoize_in
+
 
 # FIXME caches: fix up queues
 
@@ -49,12 +53,67 @@ class EvaluationMapper(EvaluationMapperBase):
 
     # {{{ map_XXX
 
+    def _map_minmax(self, func, inherited, expr):
+        ev_children = [self.rec(ch) for ch in expr.children]
+        from functools import reduce, partial
+        if any(isinstance(ch, cl.array.Array) for ch in ev_children):
+            return reduce(partial(func, queue=self.queue), ev_children)
+        else:
+            return inherited(expr)
+
+    def map_max(self, expr):
+        return self._map_minmax(
+                cl.array.maximum,
+                super(EvaluationMapper, self).map_max,
+                expr)
+
+    def map_min(self, expr):
+        return self._map_minmax(
+                cl.array.minimum,
+                super(EvaluationMapper, self).map_min,
+                expr)
+
     def map_node_sum(self, expr):
         return cl.array.sum(self.rec(expr.operand)).get()[()]
 
     def map_node_max(self, expr):
         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():
+            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)
+            knl = lp.tag_inames(knl, "el:g.0,idof:l.0")
+            return knl
+
+        discr = self.bound_expr.get_discretization(expr.where)
+
+        operand = self.rec(expr.operand)
+
+        assert operand.shape == (discr.nnodes,)
+
+        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
+
+    def map_elementwise_sum(self, expr):
+        return self._map_elementwise_reduction("sum", expr)
+
+    def map_elementwise_min(self, expr):
+        return self._map_elementwise_reduction("min", expr)
+
+    def map_elementwise_max(self, expr):
+        return self._map_elementwise_reduction("max", expr)
+
     def map_ones(self, expr):
         discr = self.bound_expr.get_discretization(expr.where)
 
@@ -73,9 +132,11 @@ class EvaluationMapper(EvaluationMapperBase):
     def map_num_reference_derivative(self, expr):
         discr = self.bound_expr.get_discretization(expr.where)
 
+        from pytools import flatten
+        ref_axes = flatten([axis] * mult for axis, mult in expr.ref_axes)
         return discr.num_reference_derivative(
                 self.queue,
-                expr.ref_axes, self.rec(expr.operand)) \
+                ref_axes, self.rec(expr.operand)) \
                         .with_queue(self.queue)
 
     def map_q_weight(self, expr):
@@ -253,6 +314,11 @@ class BoundExpression:
         return self.caches.setdefault(name, {})
 
     def get_discretization(self, where):
+        from pytential.symbolic.primitives import _QBXSourceStage2
+        if isinstance(where, _QBXSourceStage2):
+            lpot_source = self.places[where.where]
+            return lpot_source.stage2_density_discr
+
         discr = self.places[where]
 
         from pytential.source import LayerPotentialSourceBase
diff --git a/pytential/symbolic/mappers.py b/pytential/symbolic/mappers.py
index cca6b91983eeefa4c7b3e0113db671893917961f..0dfaa4f9a40a1e1afbfd3dce2f57a0bb09af4c1d 100644
--- a/pytential/symbolic/mappers.py
+++ b/pytential/symbolic/mappers.py
@@ -62,6 +62,12 @@ class IdentityMapper(IdentityMapperBase):
 
     map_node_max = map_node_sum
 
+    def map_elementwise_sum(self, expr):
+        return type(expr)(self.rec(expr.operand), expr.where)
+
+    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)
@@ -104,6 +110,9 @@ class CombineMapper(CombineMapperBase):
 
     map_node_max = map_node_sum
     map_num_reference_derivative = map_node_sum
+    map_elementwise_sum = map_node_sum
+    map_elementwise_min = map_node_sum
+    map_elementwise_max = map_node_sum
 
     def map_int_g(self, expr):
         return self.combine(
@@ -229,6 +238,14 @@ class LocationTagger(CSECachingMapperMixin, IdentityMapper):
         else:
             return expr
 
+    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)
+
+    map_elementwise_min = map_elementwise_sum
+    map_elementwise_max = map_elementwise_sum
+
     def map_int_g(self, expr):
         source = expr.source
         target = expr.target
@@ -436,6 +453,9 @@ class QBXPreprocessor(IdentityMapper):
 # {{{ stringifier
 
 def stringify_where(where):
+    if isinstance(where, prim._QBXSourceStage2):
+        return "stage2(%s)" % stringify_where(where.where)
+
     if where is None:
         return "?"
     elif where is prim.DEFAULT_SOURCE:
@@ -465,8 +485,20 @@ class StringifyMapper(BaseStringifyMapper):
                 for name_expr in six.itervalues(expr.extra_vars)),
                 set())
 
-    def map_node_sum(self, expr, enclosing_prec):
-        return "NodeSum(%s)" % self.rec(expr.operand, PREC_NONE)
+    def map_elementwise_sum(self, expr, enclosing_prec):
+        return "ElwiseSum.%s(%s)" % (
+                stringify_where(expr.where),
+                self.rec(expr.operand, PREC_NONE))
+
+    def map_elementwise_min(self, expr, enclosing_prec):
+        return "ElwiseMin.%s(%s)" % (
+                stringify_where(expr.where),
+                self.rec(expr.operand, PREC_NONE))
+
+    def map_elementwise_max(self, expr, enclosing_prec):
+        return "ElwiseMax.%s(%s)" % (
+                stringify_where(expr.where),
+                self.rec(expr.operand, PREC_NONE))
 
     def map_node_max(self, expr, enclosing_prec):
         return "NodeMax(%s)" % self.rec(expr.operand, PREC_NONE)
@@ -476,8 +508,14 @@ class StringifyMapper(BaseStringifyMapper):
                 stringify_where(expr.where))
 
     def map_num_reference_derivative(self, expr, enclosing_prec):
-        result = "d/dr%s.%s %s" % (
-                ",".join(str(ax) for ax in expr.ref_axes),
+        diff_op = " ".join(
+                "d/dr%d" % axis
+                if mult == 1 else
+                "d/dr%d^%d" % (axis, mult)
+                for axis, mult in expr.ref_axes)
+
+        result = "%s.%s %s" % (
+                diff_op,
                 stringify_where(expr.where),
                 self.rec(expr.operand, PREC_PRODUCT),
                 )
diff --git a/pytential/symbolic/primitives.py b/pytential/symbolic/primitives.py
index d981fdd2fe9e13c31093a5f67116bb6881beeba4..ecbdc89926887740e7ccd29554fa6eec491c7bc9 100644
--- a/pytential/symbolic/primitives.py
+++ b/pytential/symbolic/primitives.py
@@ -121,6 +121,10 @@ Discretization properties
 .. autofunction:: area_element
 .. autofunction:: sqrt_jac_q_weight
 .. autofunction:: normal
+.. autofunction:: mean_curvature
+.. autofunction:: first_fundamental_form
+.. autofunction:: second_fundamental_form
+.. autofunction:: shape_operator
 
 Elementary numerics
 ^^^^^^^^^^^^^^^^^^^
@@ -128,6 +132,8 @@ Elementary numerics
 .. autoclass:: NumReferenceDerivative
 .. autoclass:: NodeSum
 .. autoclass:: NodeMax
+.. autoclass:: ElementwiseSum
+.. autoclass:: ElementwiseMax
 .. autofunction:: integral
 .. autoclass:: Ones
 .. autofunction:: ones_vec
@@ -185,6 +191,8 @@ Pretty-printing expressions
 """
 
 
+# {{{ 'where' specifiers
+
 class DEFAULT_SOURCE:  # noqa
     pass
 
@@ -193,6 +201,37 @@ class DEFAULT_TARGET:  # noqa
     pass
 
 
+class _QBXSourceStage2(object):
+    """A symbolic 'where' specifier for the
+    :attr:`pytential.qbx.QBXLayerPotentialSource.stage2_density_discr`
+    of the layer potential source identified by :attr:`where`.
+
+    .. attribute:: where
+
+        An identifier of a layer potential source, as used in
+        :func:`pytential.bind`.
+
+    .. note::
+
+        This is not documented functionality and only intended for
+        internal use.
+    """
+
+    def __init__(self, where):
+        self.where = where
+
+    def __hash__(self):
+        return hash((type(self), self.where))
+
+    def __eq__(self, other):
+        return type(self) is type(other) and self.where == other.where
+
+    def __ne__(self, other):
+        return not self.__eq__(other)
+
+# }}}
+
+
 class cse_scope(cse_scope_base):  # noqa
     DISCRETIZATION = "pytential_discretization"
 
@@ -339,18 +378,46 @@ class NumReferenceDerivative(DiscretizationProperty):
     reference coordinates.
     """
 
+    def __new__(cls, ref_axes, operand, where=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 componentwise(make_op, operand)
+        else:
+            return DiscretizationProperty.__new__(cls)
+
     def __init__(self, ref_axes, operand, where=None):
         """
-        :arg ref_axes: a :class:`frozenset` of indices of
-            reference coordinates along which derivatives
-            will be taken.
+        :arg ref_axes: a :class:`tuple` of tuples indicating indices of
+            coordinate axes of the reference element to the number of derivatives
+            which will be taken.  For example, the value ``((0, 2), (1, 1))``
+            indicates that Each axis must occur at most once. The tuple must be
+            sorted by the axis index.
+
+            May also be a singile integer *i*, which is viewed as equivalent
+            to ``((i, 1),)``.
         :arg where: |where-blurb|
         """
 
-        if not isinstance(ref_axes, frozenset):
-            raise ValueError("ref_axes must be a frozenset")
+        if isinstance(ref_axes, int):
+            ref_axes = ((ref_axes, 1),)
+
+        if not isinstance(ref_axes, tuple):
+            raise ValueError("ref_axes must be a tuple")
+
+        if tuple(sorted(ref_axes)) != ref_axes:
+            raise ValueError("ref_axes must be sorted")
+
+        if len(dict(ref_axes)) != len(ref_axes):
+            raise ValueError("ref_axes must not contain an axis more than once")
 
         self.ref_axes = ref_axes
+
         self.operand = operand
         DiscretizationProperty.__init__(self, where)
 
@@ -369,10 +436,7 @@ 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(
-                frozenset([j]),
-                func_component,
-                where)
+            jac[i, j] = NumReferenceDerivative(j, func_component, where)
 
     return jac
 
@@ -382,9 +446,11 @@ def parametrization_derivative_matrix(ambient_dim, dim, where=None):
     reference-to-global parametrization.
     """
 
-    return reference_jacobian(
-            [NodeCoordinateComponent(i, where) for i in range(ambient_dim)],
-            ambient_dim, dim)
+    return cse(
+            reference_jacobian(
+                [NodeCoordinateComponent(i, where) for i in range(ambient_dim)],
+                ambient_dim, dim, where=where),
+            "pd_matrix", cse_scope.DISCRETIZATION)
 
 
 def parametrization_derivative(ambient_dim, dim, where=None):
@@ -451,9 +517,7 @@ def mean_curvature(ambient_dim, dim=None, where=None):
         raise NotImplementedError(
                 "only know how to calculate curvature for a curve in 2D")
 
-    xp, yp = cse(
-            parametrization_derivative_matrix(ambient_dim, dim, where),
-            "pd_matrix", cse_scope.DISCRETIZATION)
+    xp, yp = parametrization_derivative_matrix(ambient_dim, dim, where)
 
     xpp, ypp = cse(
             reference_jacobian([xp[0], yp[0]], ambient_dim, dim, where),
@@ -461,23 +525,222 @@ def mean_curvature(ambient_dim, dim=None, where=None):
 
     return (xp[0]*ypp[0] - yp[0]*xpp[0]) / (xp[0]**2 + yp[0]**2)**(3/2)
 
-# FIXME: make sense of this in the context of GA
-# def xyz_to_local_matrix(dim, where=None):
-#     """First two rows are tangents."""
-#     result = np.zeros((dim, dim), dtype=np.object)
-#
-#     for i in range(dim-1):
-#         result[i] = make_tangent(i, dim, where)
-#     result[-1] = make_normal(dim, where)
-#
-#     return result
+
+def first_fundamental_form(ambient_dim, dim=None, where=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)
+
+    return cse(
+            np.dot(pd_mat.T, pd_mat),
+            "fundform1")
+
+
+def second_fundamental_form(ambient_dim, dim=None, where=None):
+    """Compute the second fundamental form of a surface. This is in reference
+    to the reference-to-global mapping in use for each element.
+
+    .. note::
+
+        Some references assume that the second fundamental form is computed
+        with respect to an orthonormal basis, which this is not.
+    """
+    if dim is None:
+        dim = ambient_dim - 1
+
+    if ambient_dim != 3 and dim != 2:
+        raise NotImplementedError("only available for surfaces in 3D")
+
+    r = nodes(ambient_dim, where=where).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)
+    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()
+
+    ff2_l = cse(np.dot(ruu, nrml), "fundform2_L")
+    ff2_m = cse(np.dot(ruv, nrml), "fundform2_M")
+    ff2_n = cse(np.dot(rvv, nrml), "fundform2_N")
+
+    result = np.zeros((2, 2), dtype=object)
+    result[0, 0] = ff2_l
+    result[0, 1] = result[1, 0] = ff2_m
+    result[1, 1] = ff2_n
+
+    return result
+
+
+def shape_operator(ambient_dim, dim=None, where=None):
+    if dim is None:
+        dim = ambient_dim - 1
+
+    if ambient_dim != 3 and dim != 2:
+        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)
+
+    result = np.zeros((2, 2), dtype=object)
+    result[0, 0] = e*G-f*F
+    result[0, 1] = f*G-g*F
+    result[1, 0] = f*E-e*F
+    result[1, 1] = g*E-f*F
+
+    return cse(
+            1/(E*G-F*F)*result,
+            "shape_operator")
+
+
+def _panel_size(ambient_dim, dim=None, where=None):
+    # A broken quasi-1D approximation of 1D element size. Do not use.
+
+    if dim is None:
+        dim = ambient_dim - 1
+
+    return ElementwiseSum(
+            area_element(ambient_dim=ambient_dim, dim=dim)
+            * QWeight())**(1/dim)
+
+
+def _small_mat_inverse(mat):
+    m, n = mat.shape
+    if m != n:
+        raise ValueError("inverses only make sense for square matrices")
+
+    if m == 1:
+        return make_obj_array([1/mat[0, 0]])
+    elif m == 2:
+        (a, b), (c, d) = mat
+        return 1/(a*d-b*c) * make_obj_array([
+            [d, -b],
+            [-c, a],
+            ])
+    else:
+        raise NotImplementedError(
+                "inverse formula for %dx%d matrices" % (m, n))
+
+
+def _small_mat_eigenvalues(mat):
+    m, n = mat.shape
+    if m != n:
+        raise ValueError("eigenvalues only make sense for square matrices")
+
+    if m == 1:
+        return make_obj_array([mat[0, 0]])
+    elif m == 2:
+        (a, b), (c, d) = mat
+        return make_obj_array([
+                -(sqrt(d**2-2*a*d+4*b*c+a**2)-d-a)/2,
+                 (sqrt(d**2-2*a*d+4*b*c+a**2)+d+a)/2
+                ])
+    else:
+        raise NotImplementedError(
+                "eigenvalue formula for %dx%d matrices" % (m, n))
+
+
+def _simplex_mapping_max_stretch_factor(ambient_dim, dim=None, where=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 dim is None:
+        dim = ambient_dim - 1
+
+    # The 'technique' here is ad-hoc, but I'm fairly confident it's better than
+    # what we had. The idea is that singular values of the mapping Jacobian
+    # yield "stretch factors" of the mapping Why? because it maps a right
+    # singular vector $`v_1`$ (of unit length) to $`\sigma_1 u_1`$, where
+    # $`u_1`$ is the corresponding left singular vector (also of unit length).
+    # And so the biggest one tells us about the direction with the 'biggest'
+    # stretching, where 'stretching' (*2 to remove bi-unit reference element)
+    # reflects available quadrature resolution in that direction.
+
+    pder_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
+
+    # The above procedure works well only when the 'reference' end of the
+    # mapping is in equilateral coordinates.
+    from modepy.tools import EQUILATERAL_TO_UNIT_MAP
+    equi_to_unit = EQUILATERAL_TO_UNIT_MAP[dim].a
+
+    # This is the Jacobian of the (equilateral reference element) -> (global) map.
+    equi_pder_mat = cse(
+            np.dot(pder_mat, equi_to_unit),
+            "equilateral_pder_mat")
+
+    # Compute eigenvalues of J^T to compute SVD.
+    equi_pder_mat_jtj = cse(
+            np.dot(equi_pder_mat.T, equi_pder_mat),
+            "pd_mat_jtj")
+
+    stretch_factors = [
+            cse(sqrt(s), "mapping_singval_%d" % i)
+            for i, s in enumerate(
+                _small_mat_eigenvalues(
+                    # Multiply by 4 to compensate for equilateral reference
+                    # elements of side length 2. (J^T J contains two factors of
+                    # two.)
+                    4 * equi_pder_mat_jtj))]
+
+    from pymbolic.primitives import Max
+    result = Max(tuple(stretch_factors))
+
+    if with_elementwise_max:
+        result = ElementwiseMax(result, where=where)
+
+    return cse(result, "mapping_max_stretch", cse_scope.DISCRETIZATION)
+
+
+def _max_curvature(ambient_dim, dim=None, where=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))
+    elif ambient_dim == 3:
+        shape_op = shape_operator(ambient_dim, dim, where=where)
+
+        abs_principal_curvatures = [
+                abs(x) for x in _small_mat_eigenvalues(shape_op)]
+        from pymbolic.primitives import Max
+        return cse(Max(tuple(abs_principal_curvatures)))
+    else:
+        raise NotImplementedError("curvature criterion not implemented in %d "
+                "dimensions" % ambient_dim)
+
+
+def _scaled_max_curvature(ambient_dim, dim=None, where=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,
+                    with_elementwise_max=False)
 
 # }}}
 
 
 # {{{ operators
 
-class NodalOperation(Expression):
+class SingleScalarOperandExpression(Expression):
     def __new__(cls, operand):
         # If the constructor is handed a multivector object, return an
         # object array of the operator applied to each of the
@@ -498,13 +761,13 @@ class NodalOperation(Expression):
         return (self.operand,)
 
 
-class NodeSum(NodalOperation):
+class NodeSum(SingleScalarOperandExpression):
     """Implements a global sum over all discretization nodes."""
 
     mapper_method = "map_node_sum"
 
 
-class NodeMax(NodalOperation):
+class NodeMax(SingleScalarOperandExpression):
     """Implements a global maximum over all discretization nodes."""
 
     mapper_method = "map_node_max"
@@ -519,6 +782,52 @@ def integral(ambient_dim, dim, operand, where=None):
             * operand)
 
 
+class SingleScalarOperandExpressionWithWhere(Expression):
+    def __new__(cls, operand, where=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 componentwise(make_op, operand)
+        else:
+            return Expression.__new__(cls)
+
+    def __init__(self, operand, where=None):
+        self.operand = operand
+        self.where = where
+
+    def __getinitargs__(self):
+        return (self.operand, self.where)
+
+
+class ElementwiseSum(SingleScalarOperandExpressionWithWhere):
+    """Returns a vector of DOFs with all entries on each element set
+    to the sum of DOFs on that element.
+    """
+
+    mapper_method = "map_elementwise_sum"
+
+
+class ElementwiseMin(SingleScalarOperandExpressionWithWhere):
+    """Returns a vector of DOFs with all entries on each element set
+    to the minimum of DOFs on that element.
+    """
+
+    mapper_method = "map_elementwise_min"
+
+
+class ElementwiseMax(SingleScalarOperandExpressionWithWhere):
+    """Returns a vector of DOFs with all entries on each element set
+    to the maximum of DOFs on that element.
+    """
+
+    mapper_method = "map_elementwise_max"
+
+
 class Ones(Expression):
     """A DOF-vector that is constant *one* on the whole
     discretization.
@@ -1019,11 +1328,14 @@ def Dp(kernel, *args, **kwargs):  # noqa
 # {{{ conventional vector calculus
 
 def tangential_onb(ambient_dim, dim=None, where=None):
+    """Return a matrix of shape ``(ambient_dim, dim)`` with orthogonal columns
+    spanning the tangential space of the surface of *where*.
+    """
+
     if dim is None:
         dim = ambient_dim - 1
 
-    pd_mat = cse(parametrization_derivative_matrix(ambient_dim, dim, where),
-            "pd_matrix", cse_scope.DISCRETIZATION)
+    pd_mat = parametrization_derivative_matrix(ambient_dim, dim, where)
 
     # {{{ Gram-Schmidt
 
diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py
index 6fbd78d1082ed668eaa329afc23792dd210d420b..e7a4eaa6baa29df83e22329e4c585b7525b00b80 100644
--- a/test/test_global_qbx.py
+++ b/test/test_global_qbx.py
@@ -122,8 +122,9 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
     ext_centers = get_centers_on_side(lpot_source, +1)
     ext_centers = np.array([axis.get(queue) for axis in ext_centers])
     expansion_radii = lpot_source._expansion_radii("nsources").get(queue)
-    panel_sizes = lpot_source._panel_sizes("npanels").get(queue)
-    fine_panel_sizes = lpot_source._fine_panel_sizes("npanels").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)
 
     # {{{ check if satisfying criteria
 
@@ -155,7 +156,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
                 (dist, rad, centers_panel.element_nr, sources_panel.element_nr)
 
     def check_sufficient_quadrature_resolution(centers_panel, sources_panel):
-        h = fine_panel_sizes[sources_panel.element_nr]
+        dz_radius = source_danger_zone_radii[sources_panel.element_nr]
 
         my_int_centers = int_centers[:, centers_panel.discr_slice]
         my_ext_centers = ext_centers[:, centers_panel.discr_slice]
@@ -174,12 +175,12 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
         # Criterion:
         # The quadrature contribution from each panel is as accurate
         # as from the center's own source panel.
-        assert dist >= h / 4, \
-                (dist, h, centers_panel.element_nr, sources_panel.element_nr)
+        assert dist >= dz_radius, \
+                (dist, dz_radius, centers_panel.element_nr, sources_panel.element_nr)
 
-    def check_panel_size_to_helmholtz_k_ratio(panel):
+    def check_quad_res_to_helmholtz_k_ratio(panel):
         # Check wavenumber to panel size ratio.
-        assert panel_sizes[panel.element_nr] * helmholtz_k <= 5
+        assert quad_res[panel.element_nr] * helmholtz_k <= 5
 
     for i, panel_1 in enumerate(iter_elements(lpot_source.density_discr)):
         for panel_2 in iter_elements(lpot_source.density_discr):
@@ -187,7 +188,7 @@ def run_source_refinement_test(ctx_getter, mesh, order, helmholtz_k=None):
         for panel_2 in iter_elements(lpot_source.quad_stage2_density_discr):
             check_sufficient_quadrature_resolution(panel_1, panel_2)
         if helmholtz_k is not None:
-            check_panel_size_to_helmholtz_k_ratio(panel_1)
+            check_quad_res_to_helmholtz_k_ratio(panel_1)
 
     # }}}
 
@@ -216,10 +217,17 @@ def test_source_refinement_3d(ctx_getter, surface_name, surface_f, order):
 
 
 @pytest.mark.parametrize(("curve_name", "curve_f", "nelements"), [
-    ("20-to-1 ellipse", partial(ellipse, 20), 100),
+    # TODO: This used to pass for a 20-to-1 ellipse (with different center
+    # placement).  It still produces valid output right now, but the current
+    # test is not smart enough to recognize it. It might be useful to fix that.
+    #
+    # See discussion at
+    # https://gitlab.tiker.net/inducer/pytential/merge_requests/95#note_24003
+    ("18-to-1 ellipse", partial(ellipse, 18), 100),
     ("horseshoe", horseshoe, 64),
     ])
-def test_target_association(ctx_getter, curve_name, curve_f, nelements):
+def test_target_association(ctx_getter, curve_name, curve_f, nelements,
+        visualize=False):
     cl_ctx = ctx_getter()
     queue = cl.CommandQueue(cl_ctx)
 
@@ -315,6 +323,35 @@ def test_target_association(ctx_getter, curve_name, curve_f, nelements):
     int_targets = np.array([axis.get(queue) for axis in int_targets.nodes()])
     ext_targets = np.array([axis.get(queue) for axis in ext_targets.nodes()])
 
+    def visualize_curve_and_assoc():
+        import matplotlib.pyplot as plt
+        from meshmode.mesh.visualization import draw_curve
+
+        draw_curve(lpot_source.density_discr.mesh)
+
+        targets = int_targets
+        tgt_slice = surf_int_slice
+
+        plt.plot(centers[0], centers[1], "+", color="orange")
+        ax = plt.gca()
+
+        for tx, ty, tcenter in zip(
+                targets[0, tgt_slice],
+                targets[1, tgt_slice],
+                target_assoc.target_to_center[tgt_slice]):
+            if tcenter >= 0:
+                ax.add_artist(
+                        plt.Line2D(
+                            (tx, centers[0, tcenter]),
+                            (ty, centers[1, tcenter]),
+                            ))
+
+        ax.set_aspect("equal")
+        plt.show()
+
+    if visualize:
+        visualize_curve_and_assoc()
+
     # Checks that the sources match with their own centers.
     def check_on_surface_targets(nsources, true_side, target_to_center,
                                  target_to_side_result):
diff --git a/test/test_layer_pot.py b/test/test_layer_pot.py
index b0766375186a5ad37f9fcd5e3ef4a04fa907c964..f81a5e550d06e879dd2c9309813d97381bef97b6 100644
--- a/test/test_layer_pot.py
+++ b/test/test_layer_pot.py
@@ -404,7 +404,7 @@ def test_3d_jump_relations(ctx_factory, relation, visualize=False):
     from pytools.convergence import EOCRecorder
     eoc_rec = EOCRecorder()
 
-    for nel_factor in [6, 8, 12]:
+    for nel_factor in [6, 10, 14]:
         from meshmode.mesh.generation import generate_torus
         mesh = generate_torus(
                 5, 2, order=target_order,
diff --git a/test/test_layer_pot_eigenvalues.py b/test/test_layer_pot_eigenvalues.py
index e60b72bfc4108cf229114ff5d9c3563afa879d23..b8c9c9fd0a2ed2f13408069688bb7746724eff9e 100644
--- a/test/test_layer_pot_eigenvalues.py
+++ b/test/test_layer_pot_eigenvalues.py
@@ -100,7 +100,7 @@ def test_ellipse_eigenvalues(ctx_getter, ellipse_aspect, mode_nr, qbx_order,
                 np.linspace(0, 1, nelements+1),
                 target_order)
 
-        fmm_order = 10
+        fmm_order = 12
         if force_direct:
             fmm_order = False
 
diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py
index c019f9309f2d23b5fe910888e1229cc850dfb51f..942b6ed282ae1eb250c5d6bea645bc834dda84eb 100644
--- a/test/test_layer_pot_identity.py
+++ b/test/test_layer_pot_identity.py
@@ -82,7 +82,7 @@ class StarfishGeometry(object):
 
     dim = 2
 
-    resolutions = [30, 50, 70]
+    resolutions = [30, 50, 70, 90]
 
     def get_mesh(self, nelements, target_order):
         return make_curve_mesh(
diff --git a/test/test_scalar_int_eq.py b/test/test_scalar_int_eq.py
index f10f471e7c307e752e8a0fe1715fad7c958c132d..0d9fca2871c95f7070876e560f9cf94bafaccbc5 100644
--- a/test/test_scalar_int_eq.py
+++ b/test/test_scalar_int_eq.py
@@ -83,7 +83,7 @@ class IntEqTestCase:
 
 
 class CurveIntEqTestCase(IntEqTestCase):
-    resolutions = [30, 40, 50]
+    resolutions = [40, 50, 60]
 
     def get_mesh(self, resolution, target_order):
         return make_curve_mesh(