diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py
index dc4a3b3c6..29c31c18e 100644
--- a/meshmode/mesh/visualization.py
+++ b/meshmode/mesh/visualization.py
@@ -27,6 +27,7 @@
 .. autofunction:: draw_curve
 .. autofunction:: write_vertex_vtk_file
 .. autofunction:: mesh_to_tikz
+.. autofunction:: vtk_visualize_mesh
 """
 
 
@@ -299,4 +300,23 @@ def mesh_to_tikz(mesh):
 
 # }}}
 
+
+# {{{ visualize_mesh
+
+def vtk_visualize_mesh(actx, mesh, filename, vtk_high_order=True):
+    order = max(mgrp.order for mgrp in mesh.groups)
+
+    from meshmode.discretization.poly_element import \
+            PolynomialWarpAndBlendGroupFactory
+    from meshmode.discretization import Discretization
+    discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+
+    from meshmode.discretization.visualization import make_visualizer
+    vis = make_visualizer(actx, discr, order, force_equidistant=vtk_high_order)
+
+    vis.write_vtk_file(filename, [], use_high_order=vtk_high_order)
+
+# }}}
+
+
 # vim: foldmethod=marker
diff --git a/test/test_refinement.py b/test/test_refinement.py
index c88f7a148..1835edb3a 100644
--- a/test/test_refinement.py
+++ b/test/test_refinement.py
@@ -36,6 +36,7 @@
         generate_icosahedron, generate_box_mesh, make_curve_mesh, ellipse)
 from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry
 from meshmode.mesh.refinement import Refiner, RefinerWithoutAdjacency
+import meshmode.mesh.generation as mgen
 
 from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup
 from meshmode.discretization.poly_element import (
@@ -224,8 +225,7 @@ def test_refinement_connection(
             mesh = get_blob_mesh(mesh_par, mesh_order)
             h = float(mesh_par)
         elif mesh_name == "warp":
-            from meshmode.mesh.generation import generate_warped_rect_mesh
-            mesh = generate_warped_rect_mesh(dim, order=mesh_order, n=mesh_par,
+            mesh = mgen.generate_warped_rect_mesh(dim, order=mesh_order, n=mesh_par,
                     group_cls=group_cls)
             h = 1/mesh_par
         else:
@@ -321,8 +321,7 @@ def test_uniform_refinement(group_cls, with_adjacency):
 
 @pytest.mark.parametrize("refinement_rounds", [0, 1, 2])
 def test_conformity_of_uniform_mesh(refinement_rounds):
-    from meshmode.mesh.generation import generate_icosphere
-    mesh = generate_icosphere(r=1.0, order=4,
+    mesh = mgen.generate_icosphere(r=1.0, order=4,
             uniform_refinement_rounds=refinement_rounds)
 
     assert mesh.is_conforming
@@ -331,6 +330,30 @@ def test_conformity_of_uniform_mesh(refinement_rounds):
     assert is_boundary_tag_empty(mesh, BTAG_ALL)
 
 
+@pytest.mark.parametrize("mesh_name", ["torus", "icosphere"])
+def test_refine_surfaces(actx_factory, mesh_name, visualize=False):
+    if mesh_name == "torus":
+        mesh = mgen.generate_torus(10, 1, 40, 4, order=4)
+    elif mesh_name == "icosphere":
+        mesh = mgen.generate_icosphere(1, order=4)
+    else:
+        raise ValueError(f"invalid mesh name '{mesh_name}'")
+
+    if visualize:
+        actx = actx_factory()
+        from meshmode.mesh.visualization import vtk_visualize_mesh
+        vtk_visualize_mesh(actx, mesh, "surface.vtu")
+
+    # check for absence of node-vertex consistency error
+    from meshmode.mesh.refinement import refine_uniformly
+    refined_mesh = refine_uniformly(mesh, 1)
+
+    if visualize:
+        actx = actx_factory()
+        from meshmode.mesh.visualization import vtk_visualize_mesh
+        vtk_visualize_mesh(actx, refined_mesh, "surface-refined.vtu")
+
+
 if __name__ == "__main__":
     import sys
     if len(sys.argv) > 1: