diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 0b0c55985..b1fb866ed 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -70,7 +70,7 @@ jobs:
                 grep -v symengine .test-conda-env-py3.yml > .test-conda-env.yml
                 CONDA_ENVIRONMENT=.test-conda-env.yml
 
-                export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:--k-slowtest}
+                export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-k 'not slowtest'"}
 
                 curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh
                 . ./build-and-test-py-project-within-miniconda.sh
@@ -88,7 +88,7 @@ jobs:
                 echo "- compilers" >> .test-conda-env.yml
                 echo "- llvm-openmp" >> .test-conda-env.yml
                 CONDA_ENVIRONMENT=.test-conda-env.yml
-                export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:--k-slowtest}
+                export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-k 'not slowtest'"}
 
                 curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh
                 . ./build-and-test-py-project-within-miniconda.sh
@@ -102,7 +102,7 @@ jobs:
             run: |
                 CONDA_ENVIRONMENT=.test-conda-env-py3.yml
 
-                export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:--k-slowtest}
+                export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-k 'not slowtest'"}
 
                 curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh
                 . ./build-and-test-py-project-within-miniconda.sh
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b78c640e5..ad92ccd24 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,6 +1,6 @@
 # Environment variables
 #
-# * PYTEST_ADDOPTS is used to filter test runs. The default value is "-k-slowtest",
+# * PYTEST_ADDOPTS is used to filter test runs. The default value is "-k 'not slowtest'",
 #   which skips the slow running tests.
 # * SKIP_EXAMPLES, if non-empty, can be used to skip the examples job.
 
@@ -8,7 +8,7 @@ Python 3 POCL:
   script:
   - export PY_EXE=python3
   - export PYOPENCL_TEST=portable:pthread
-  - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:--k-slowtest}
+  - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-k 'not slowtest'"}
   - export EXTRA_INSTALL="Cython pybind11 numpy scipy mako"
   - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh
   - ". ./build-and-test-py-project.sh"
@@ -27,7 +27,7 @@ Python 3 Intel:
   - export PY_EXE=python3
   - source /opt/enable-intel-cl.sh
   - export PYOPENCL_TEST="intel(r):pu"
-  - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:--k-slowtest}
+  - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-k 'not slowtest'"}
   - export EXTRA_INSTALL="Cython pybind11 numpy scipy mako"
   - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project.sh
   - ". ./build-and-test-py-project.sh"
@@ -60,7 +60,7 @@ Python 3 Conda:
   script:
   - export SUMPY_FORCE_SYMBOLIC_BACKEND=symengine
   - export CONDA_ENVIRONMENT=.test-conda-env-py3.yml
-  - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:--k-slowtest}
+  - export PYTEST_ADDOPTS=${PYTEST_ADDOPTS:-"-k 'not slowtest'"}
 
   - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/main/build-and-test-py-project-within-miniconda.sh
   - ". ./build-and-test-py-project-within-miniconda.sh"
diff --git a/experiments/maxwell_sphere.py b/experiments/maxwell_sphere.py
index b2ce1940e..d1101ff7d 100644
--- a/experiments/maxwell_sphere.py
+++ b/experiments/maxwell_sphere.py
@@ -26,17 +26,8 @@
 def main():
     # cl.array.to_device(queue, numpy_array)
     from meshmode.mesh.io import generate_gmsh, FileSource
-    from meshmode.mesh.generation import generate_icosphere
-    from meshmode.mesh.refinement import Refiner
-    mesh = generate_icosphere(1,target_order)
-
-    refinement_increment = 1
-    refiner = Refiner(mesh)
-    for i in range(refinement_increment):
-        flags = np.ones(mesh.nelements, dtype=bool)
-        refiner.refine(flags)
-        mesh = refiner.get_current_mesh()
-
+    from meshmode.mesh.generation import generate_sphere
+    mesh = generate_sphere(1, target_order, uniform_refinement_rounds=1)
 
     from meshmode.mesh.processing import perform_flips
     # Flip elements--gmsh generates inside-out geometry.
@@ -135,7 +126,7 @@ def green3m(x,y,z,source,strength,k):
         # magnetic field corresponding to dyadic green's function
         # due to monochromatic electric dipole located at "source".
         # "strength" is the the intensity of the dipole.
-        #  H = curl((I + Hess)(exp(ikr)/r) dot (strength)) = 
+        #  H = curl((I + Hess)(exp(ikr)/r) dot (strength)) =
         #  strength \cross \grad (exp(ikr)/r)
         #
             dx = x - source[0]
@@ -165,7 +156,7 @@ def dipole3e(x,y,z,source,strength,k):
 #            print(hvec)
 #            print(strength)
             return evec,hvec
-            
+
         def dipole3m(x,y,z,source,strength,k):
         #
         #  evalaute electric and magnetic field due
@@ -175,7 +166,7 @@ def dipole3m(x,y,z,source,strength,k):
             hvec = green3e(x,y,z,source,strength,k)
             hvec = -hvec*1j*k
             return evec,hvec
-            
+
 
         def dipole3eall(x,y,z,sources,strengths,k):
             ns = len(strengths)
@@ -194,7 +185,7 @@ def dipole3eall(x,y,z,sources,strengths,k):
 #        source[1] =-0.03
 #        source[2] = 0.02
         strength = np.ones(3)
-       
+
 #        evec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128))
 #        hvec = cl.array.to_device(queue,np.zeros((3,len(nodes[0])),dtype=np.complex128))
 
diff --git a/pytential/linalg/proxy.py b/pytential/linalg/proxy.py
index 4c8b58dff..1af52e4a0 100644
--- a/pytential/linalg/proxy.py
+++ b/pytential/linalg/proxy.py
@@ -451,8 +451,7 @@ def prg():
                 <> rqbx_int = simul_reduce(max, i, sqrt(simul_reduce(sum, idim, \
                         (proxy_centers[idim, irange] -
                          center_int[idim, srcindices[i + ioffset]]) ** 2)) + \
-                         expansion_radii[srcindices[i + ioffset]]) \
-                         {dup=idim}
+                         expansion_radii[srcindices[i + ioffset]])
                 <> rqbx_ext = simul_reduce(max, i, sqrt(simul_reduce(sum, idim, \
                         (proxy_centers[idim, irange] -
                          center_ext[idim, srcindices[i + ioffset]]) ** 2)) + \
diff --git a/pytential/solve.py b/pytential/solve.py
index 979a7550b..b2e819665 100644
--- a/pytential/solve.py
+++ b/pytential/solve.py
@@ -31,27 +31,32 @@
 .. autoclass:: ResidualPrinter
 """
 
+from functools import partial
 
 import numpy as np
-from numbers import Number
-import pyopencl as cl
-import pyopencl.array  # noqa
-from meshmode.dof_array import obj_or_dof_array_vectorize_n_args, DOFArray
-from pytools.obj_array import obj_array_vectorize_n_args
 
 
-def structured_vdot(x, y):
-    # vdot() implementation that is aware of scalars and host or
-    # PyOpenCL arrays. It also recurses down nested object arrays.
+def structured_vdot(x, y, array_context=None):
+    """vdot() implementation that is aware of scalars and host or
+    PyOpenCL arrays. It also recurses down nested object arrays.
+    """
+
+    if not type(x) == type(y):
+        raise TypeError("'structured_vdot' entries have different types: "
+                f"{type(x).__name__} and {type(y).__name__}")
+
+    from numbers import Number
     if (isinstance(x, Number)
             or (isinstance(x, np.ndarray) and x.dtype.char != "O")):
         return np.vdot(x, y)
-    elif isinstance(x, cl.array.Array):
-        return cl.array.vdot(x, y).get()
-    elif isinstance(x, np.ndarray) and x.dtype.char == "O":
-        return np.sum(obj_array_vectorize_n_args(structured_vdot, x, y))
-    elif isinstance(x, DOFArray):
-        return sum(obj_or_dof_array_vectorize_n_args(structured_vdot, x, y))
+    else:
+        # actx.np.vdot works on PyOpenCL arrays and arbitrarily nested
+        # array containers, so this should handle all remaining cases
+        r = array_context.to_numpy(array_context.np.vdot(x, y))
+        if isinstance(r, np.ndarray) and r.shape == ():
+            r = r[()]
+
+        return r
 
 
 # {{{ gmres
@@ -262,7 +267,7 @@ def __call__(self, resid):
 # {{{ entrypoint
 
 def gmres(op, rhs, restart=None, tol=None, x0=None,
-        inner_product=structured_vdot,
+        inner_product=None,
         maxiter=None, hard_failure=None,
         no_progress_factor=None, stall_iterations=None,
         callback=None, progress=False, require_monotonicity=True):
@@ -284,6 +289,18 @@ def gmres(op, rhs, restart=None, tol=None, x0=None,
 
     :return: a :class:`GMRESResult`
     """
+    if inner_product is None:
+        from pytential.symbolic.execution import \
+                _find_array_context_from_args_in_context
+        try:
+            actx = _find_array_context_from_args_in_context({
+                "rhs": rhs, "x0": x0,
+                }, supplied_array_context=getattr(op, "array_context", None))
+        except (ValueError, TypeError):
+            actx = None
+
+        inner_product = partial(structured_vdot, array_context=actx)
+
     if callback is None:
         if progress:
             callback = ResidualPrinter(inner_product)
diff --git a/pytential/symbolic/execution.py b/pytential/symbolic/execution.py
index 0984bf0d8..5695d3987 100644
--- a/pytential/symbolic/execution.py
+++ b/pytential/symbolic/execution.py
@@ -1011,9 +1011,6 @@ def cost_per_stage(self, calibration_params, **kwargs):
         """
         array_context = _find_array_context_from_args_in_context(kwargs)
 
-        if array_context is None:
-            raise ValueError("unable to figure array context from arguments")
-
         cost_model_mapper = CostModelMapper(
             self, array_context, calibration_params, per_box=False, context=kwargs
         )
diff --git a/pytential/symbolic/matrix.py b/pytential/symbolic/matrix.py
index 65a8579a4..949bf9706 100644
--- a/pytential/symbolic/matrix.py
+++ b/pytential/symbolic/matrix.py
@@ -346,10 +346,14 @@ def map_interpolation(self, expr):
             except KeyError:
                 from meshmode.discretization.connection import \
                     flatten_chained_connection
+                from meshmode.discretization.connection.direct import \
+                    make_direct_full_resample_matrix
 
                 conn = self.places.get_connection(expr.from_dd, expr.to_dd)
                 conn = flatten_chained_connection(actx, conn)
-                mat = actx.to_numpy(conn.full_resample_matrix(actx))
+                mat = actx.to_numpy(
+                    make_direct_full_resample_matrix(actx, conn)
+                    )
 
                 # FIXME: the resample matrix is slow to compute and very big
                 # to store, so caching it may not be the best idea
diff --git a/test/extra_int_eq_data.py b/test/extra_int_eq_data.py
index d3cecdacf..c1f0893fd 100644
--- a/test/extra_int_eq_data.py
+++ b/test/extra_int_eq_data.py
@@ -362,8 +362,8 @@ class SphereTestCase(IntegralEquationTestCase):
     check_tangential_deriv = False
 
     def get_mesh(self, resolution, mesh_order):
-        from meshmode.mesh.generation import generate_icosphere
-        return generate_icosphere(1.0, mesh_order,
+        from meshmode.mesh.generation import generate_sphere
+        return generate_sphere(1.0, mesh_order,
                 uniform_refinement_rounds=resolution)
 
 
diff --git a/test/test_cost_model.py b/test/test_cost_model.py
index 1f389aeb6..1402e3bec 100644
--- a/test/test_cost_model.py
+++ b/test/test_cost_model.py
@@ -347,7 +347,7 @@ def test_timing_data_gathering(ctx_factory):
     pytest.importorskip("pyfmmlib")
 
     import pyopencl as cl
-    from arraycontext import PyOpenCLArrayContext
+    from meshmode.array_context import PyOpenCLArrayContext
     cl_ctx = ctx_factory()
     queue = cl.CommandQueue(cl_ctx,
             properties=cl.command_queue_properties.PROFILING_ENABLE)
diff --git a/test/test_global_qbx.py b/test/test_global_qbx.py
index 2bd8a41b4..65c79debe 100644
--- a/test/test_global_qbx.py
+++ b/test/test_global_qbx.py
@@ -235,7 +235,7 @@ def test_source_refinement_2d(actx_factory, curve_name, curve_f, nelements):
 
 
 @pytest.mark.parametrize(("surface_name", "surface_f", "order"), [
-    ("sphere", partial(mgen.generate_icosphere, 1), 4),
+    ("sphere", partial(mgen.generate_sphere, 1), 4),
     ("torus", partial(mgen.generate_torus, 3, 1, n_minor=10, n_major=7), 6),
     ])
 def test_source_refinement_3d(actx_factory, surface_name, surface_f, order):
diff --git a/test/test_layer_pot_eigenvalues.py b/test/test_layer_pot_eigenvalues.py
index d738a9e3c..2ec593ce2 100644
--- a/test/test_layer_pot_eigenvalues.py
+++ b/test/test_layer_pot_eigenvalues.py
@@ -279,15 +279,9 @@ def rel_err(comp, ref):
                 )
 
     for nrefinements in [0, 1]:
-        from meshmode.mesh.generation import generate_icosphere
-        mesh = generate_icosphere(1, target_order)
-        from meshmode.mesh.refinement import Refiner
-
-        refiner = Refiner(mesh)
-        for _ in range(nrefinements):
-            flags = np.ones(mesh.nelements, dtype=bool)
-            refiner.refine(flags)
-            mesh = refiner.get_current_mesh()
+        from meshmode.mesh.generation import generate_sphere
+        mesh = generate_sphere(1, target_order,
+                uniform_refinement_rounds=nrefinements)
 
         pre_density_discr = Discretization(
                 actx, mesh,
diff --git a/test/test_layer_pot_identity.py b/test/test_layer_pot_identity.py
index 21ae59b69..e30d2a5a9 100644
--- a/test/test_layer_pot_identity.py
+++ b/test/test_layer_pot_identity.py
@@ -54,17 +54,9 @@
 
 
 def get_sphere_mesh(refinement_increment, target_order):
-    from meshmode.mesh.generation import generate_icosphere
-    mesh = generate_icosphere(1, target_order)
-    from meshmode.mesh.refinement import Refiner
-
-    refiner = Refiner(mesh)
-    for _ in range(refinement_increment):
-        flags = np.ones(mesh.nelements, dtype=bool)
-        refiner.refine(flags)
-        mesh = refiner.get_current_mesh()
-
-    return mesh
+    from meshmode.mesh.generation import generate_sphere
+    return generate_sphere(1, target_order,
+            uniform_refinement_rounds=refinement_increment)
 
 
 class StarfishGeometry:
diff --git a/test/test_linalg_proxy.py b/test/test_linalg_proxy.py
index 4e9a59138..425400923 100644
--- a/test/test_linalg_proxy.py
+++ b/test/test_linalg_proxy.py
@@ -152,8 +152,8 @@ def plot_proxy_geometry(
             from meshmode.discretization.poly_element import \
                 InterpolatoryQuadratureSimplexGroupFactory
 
-            from meshmode.mesh.generation import generate_icosphere
-            ref_mesh = generate_icosphere(1, 4, uniform_refinement_rounds=1)
+            from meshmode.mesh.generation import generate_sphere
+            ref_mesh = generate_sphere(1, 4, uniform_refinement_rounds=1)
             pxycenters = np.stack(pxy.centers)
 
             for i in range(indices.nblocks):
diff --git a/test/test_maxwell.py b/test/test_maxwell.py
index 590e92ff5..50288443c 100644
--- a/test/test_maxwell.py
+++ b/test/test_maxwell.py
@@ -65,19 +65,17 @@ class SphereTestCase(MaxwellTestCase):
     gmres_tol = 1e-10
 
     def get_mesh(self, resolution, target_order):
-        from meshmode.mesh.generation import generate_icosphere
-        from meshmode.mesh.refinement import refine_uniformly
-        return refine_uniformly(
-                generate_icosphere(2, target_order),
-                resolution)
+        from meshmode.mesh.generation import generate_sphere
+        return generate_sphere(2, target_order,
+                uniform_refinement_rounds=resolution)
 
     def get_observation_mesh(self, target_order):
-        from meshmode.mesh.generation import generate_icosphere
+        from meshmode.mesh.generation import generate_sphere
 
         if self.is_interior:
-            return generate_icosphere(5, target_order)
+            return generate_sphere(5, target_order)
         else:
-            return generate_icosphere(0.5, target_order)
+            return generate_sphere(0.5, target_order)
 
     def get_source(self, actx):
         if self.is_interior:
@@ -114,12 +112,12 @@ def get_mesh(self, resolution, target_order):
         return perform_flips(mesh, np.ones(mesh.nelements))
 
     def get_observation_mesh(self, target_order):
-        from meshmode.mesh.generation import generate_icosphere
+        from meshmode.mesh.generation import generate_sphere
 
         if self.is_interior:
-            return generate_icosphere(5, target_order)
+            return generate_sphere(5, target_order)
         else:
-            return generate_icosphere(0.5, target_order)
+            return generate_sphere(0.5, target_order)
 
     def get_source(self, actx):
         if self.is_interior:
@@ -159,12 +157,12 @@ def get_mesh(self, resolution, target_order):
         return perform_flips(mesh, np.ones(mesh.nelements))
 
     def get_observation_mesh(self, target_order):
-        from meshmode.mesh.generation import generate_icosphere
+        from meshmode.mesh.generation import generate_sphere
 
         if self.is_interior:
-            return generate_icosphere(12, target_order)
+            return generate_sphere(12, target_order)
         else:
-            return generate_icosphere(0.5, target_order)
+            return generate_sphere(0.5, target_order)
 
     def get_source(self, actx):
         if self.is_interior:
diff --git a/test/test_stokes.py b/test/test_stokes.py
index 666bac010..5c08544c2 100644
--- a/test/test_stokes.py
+++ b/test/test_stokes.py
@@ -70,8 +70,8 @@ def run_exterior_stokes(actx_factory, *,
                 np.linspace(0.0, 1.0, resolution + 1),
                 target_order)
     elif ambient_dim == 3:
-        from meshmode.mesh.generation import generate_icosphere
-        mesh = generate_icosphere(radius, target_order + 1,
+        from meshmode.mesh.generation import generate_sphere
+        mesh = generate_sphere(radius, target_order + 1,
                 uniform_refinement_rounds=resolution)
     else:
         raise ValueError(f"unsupported dimension: {ambient_dim}")
diff --git a/test/test_target_specific_qbx.py b/test/test_target_specific_qbx.py
index f89758750..0664a0b6c 100644
--- a/test/test_target_specific_qbx.py
+++ b/test/test_target_specific_qbx.py
@@ -137,8 +137,8 @@ def test_target_specific_qbx(actx_factory, op, helmholtz_k, qbx_order):
     target_order = 4
     fmm_tol = 1e-3
 
-    from meshmode.mesh.generation import generate_icosphere
-    mesh = generate_icosphere(1, target_order)
+    from meshmode.mesh.generation import generate_sphere
+    mesh = generate_sphere(1, target_order)
 
     from meshmode.discretization import Discretization
     from meshmode.discretization.poly_element import \