diff --git a/pytential/qbx/__init__.py b/pytential/qbx/__init__.py
index ff4c1adece818fa0068c47030456501bd899b337..d3cd843a7e1a0b42ab49be23e02d3e6c35d302ef 100644
--- a/pytential/qbx/__init__.py
+++ b/pytential/qbx/__init__.py
@@ -375,6 +375,28 @@ class QBXLayerPotentialSource(LayerPotentialSourceBase):
 
         return conn
 
+    @property
+    @memoize_method
+    def direct_resampler(self):
+        """
+        .. warning::
+
+            This always returns a
+            :class:`~meshmode.discretization.connection.DirectDiscretizationConnect`.
+            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
+            the one returned by :attr:`resampler`.
+        """
+        from meshmode.discretization.connection import \
+                flatten_chained_connection
+
+        conn = self.resampler
+        with cl.CommandQueue(self.cl_context) as queue:
+            conn = flatten_chained_connection(queue, conn)
+
+        return conn
+
     @property
     @memoize_method
     def tree_code_container(self):
diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py
index 84cd46ce23990fb893f47f3e3c2d7672aa40129f..f0b1d82f7ae9c963341450e07bb59ba761df9967 100644
--- a/pytential/symbolic/matrix.py
+++ b/pytential/symbolic/matrix.py
@@ -198,8 +198,9 @@ class MatrixBuilder(EvaluationMapperBase):
         waa = source.weights_and_area_elements().get(queue=self.queue)
         mat[:, :] *= waa
 
-        resample_mat = (
-                source.resampler.full_resample_matrix(self.queue).get(self.queue))
+        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)
 
diff --git a/test/test_matrix.py b/test/test_matrix.py
index 1d04789dbc254df1e9b67b6fbec48717a642af05..9f0ff0b44055dcb0ff91eda7e410c0a0f9149e65 100644
--- a/test/test_matrix.py
+++ b/test/test_matrix.py
@@ -26,9 +26,8 @@ import numpy as np
 import numpy.linalg as la
 import pyopencl as cl
 import pytest
-from meshmode.mesh.generation import (  # noqa
-        ellipse, cloverleaf, starfish, drop, n_gon, qbx_peanut,
-        make_curve_mesh)
+from meshmode.mesh.generation import \
+        ellipse, NArmedStarfish, make_curve_mesh
 from pytential import bind, sym
 from functools import partial
 from sumpy.symbolic import USE_SYMENGINE
@@ -40,10 +39,12 @@ from pyopencl.tools import (  # noqa
 
 @pytest.mark.skipif(USE_SYMENGINE,
         reason="https://gitlab.tiker.net/inducer/sumpy/issues/25")
-@pytest.mark.parametrize(("k", "layer_pot_id"),
-                        [(0, 1), (0, 2),
-                         (42, 1), (42, 2)])
-def test_matrix_build(ctx_factory, k, layer_pot_id, visualize=False):
+@pytest.mark.parametrize("k", [0, 42])
+@pytest.mark.parametrize("curve_f", [
+    partial(ellipse, 3),
+    NArmedStarfish(5, 0.25)])
+@pytest.mark.parametrize("layer_pot_id", [1, 2])
+def test_matrix_build(ctx_factory, k, curve_f, layer_pot_id, visualize=False):
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx)
 
@@ -54,7 +55,6 @@ def test_matrix_build(ctx_factory, k, layer_pot_id, visualize=False):
     target_order = 7
     qbx_order = 4
     nelements = 30
-    curve_f = partial(ellipse, 3)
 
     from sumpy.kernel import LaplaceKernel, HelmholtzKernel
     if k:
@@ -106,7 +106,8 @@ def test_matrix_build(ctx_factory, k, layer_pot_id, visualize=False):
 
     if visualize:
         from sumpy.tools import build_matrix as build_matrix_via_matvec
-        mat2 = build_matrix_via_matvec(bound_op.scipy_op(queue, "u"))
+        mat2 = bound_op.scipy_op(queue, "u", dtype=mat.dtype, **knl_kwargs)
+        mat2 = build_matrix_via_matvec(mat2)
         print(la.norm((mat - mat2).real, "fro") / la.norm(mat2.real, "fro"),
               la.norm((mat - mat2).imag, "fro") / la.norm(mat2.imag, "fro"))
 
@@ -135,7 +136,7 @@ def test_matrix_build(ctx_factory, k, layer_pot_id, visualize=False):
         if is_obj_array(u_sym):
             u = make_obj_array([
                 np.random.randn(density_discr.nnodes)
-                for i in range(len(u_sym))
+                for _ in range(len(u_sym))
                 ])
         else:
             u = np.random.randn(density_discr.nnodes)