diff --git a/examples/refinement-playground.py b/examples/refinement-playground.py
index 91daa3f07145037aba91ee33df874fa4818bd20e..a29c0db8f687440f3493015cb256d30926f6b4d8 100644
--- a/examples/refinement-playground.py
+++ b/examples/refinement-playground.py
@@ -1,7 +1,5 @@
 from __future__ import division, print_function
 
-from meshmode.mesh.refinement.utils import test_nodal_adj_against_geometry
-
 from six.moves import range
 import numpy as np  # noqa
 import pyopencl as cl
@@ -10,6 +8,9 @@ import os
 import logging
 order = 1
 
+from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry
+from meshmode.mesh.refinement import Refiner
+
 
 #construct vertex vertex_index
 def remove_if_exists(name):
@@ -116,11 +117,9 @@ def refine_and_generate_chart_function(mesh, filename, function):
         for k in range(len(mesh.vertices)):
             print(mesh.vertices[k, i])
 
-    #test_nodal_adj_against_geometry(mesh);
-    from meshmode.mesh.refinement import Refiner
+    #check_nodal_adj_against_geometry(mesh);
     r = Refiner(mesh)
     #random.seed(0)
-    poss_flags = np.ones(len(mesh.groups[0].vertex_indices))
     #times = 3
     num_elements = []
     time_t = []
@@ -175,7 +174,7 @@ def refine_and_generate_chart_function(mesh, filename, function):
     flags[2] = 1
     mesh = r.refine(flags)
     '''
-    #test_nodal_adj_against_geometry(mesh)
+    #check_nodal_adj_against_geometry(mesh)
     #r.print_rays(70)
     #r.print_rays(117)
     #r.print_hanging_elements(10)
@@ -231,7 +230,6 @@ def all_refine(num_mesh, depth, fname):
     runtimes = []
     for el_fact in range(2, num_mesh+2):
         mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
-        from meshmode.mesh.refinement import Refiner
         r = Refiner(mesh)
         for time in range(depth):
             flags = np.ones(len(mesh.groups[0].vertex_indices))
@@ -243,7 +241,7 @@ def all_refine(num_mesh, depth, fname):
                 stop = timeit.default_timer()
                 nelements.append(mesh.nelements)
                 runtimes.append(stop-start)
-        test_nodal_adj_against_geometry(mesh)
+        check_nodal_adj_against_geometry(mesh)
     import matplotlib.pyplot as pt
     pt.plot(nelements, runtimes, "o")
     pt.savefig(fname)
@@ -261,7 +259,6 @@ def uniform_refine(num_mesh, fract, depth, fname):
     runtimes = []
     for el_fact in range(2, num_mesh+2):
         mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),))
-        from meshmode.mesh.refinement import Refiner
         r = Refiner(mesh)
         all_els = list(range(mesh.nelements))
         for time in range(depth):
@@ -291,7 +288,7 @@ def uniform_refine(num_mesh, fract, depth, fname):
                     all_els.append(i)
             for i in range(len(flags), mesh.nelements):
                 all_els.append(i)
-            test_nodal_adj_against_geometry(mesh)
+            check_nodal_adj_against_geometry(mesh)
 
     import matplotlib.pyplot as pt
     pt.plot(nelements, runtimes, "o")
diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py
index 7f22fa5cfd155016ccab36db8aaf8fbdbaf74d5b..9a78c31d6cd2dbe5bdb2f3729815b06630c96c7f 100644
--- a/meshmode/mesh/refinement/utils.py
+++ b/meshmode/mesh/refinement/utils.py
@@ -43,7 +43,7 @@ def is_symmetric(relation, debug=False):
     return True
 
 
-def test_nodal_adj_against_geometry(mesh, tol=1e-12):
+def check_nodal_adj_against_geometry(mesh, tol=1e-12):
     def group_and_iel_to_global_iel(igrp, iel):
         return mesh.groups[igrp].element_nr_base + iel
 
diff --git a/test/test_refinement.py b/test/test_refinement.py
new file mode 100644
index 0000000000000000000000000000000000000000..aa1ab68a725c73c8618e32fa7c7a8dd518a4420b
--- /dev/null
+++ b/test/test_refinement.py
@@ -0,0 +1,135 @@
+from __future__ import division, absolute_import, print_function
+
+__copyright__ = "Copyright (C) 2014-6 Shivam Gupta, Andreas Kloeckner"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+import pytest
+
+import numpy as np
+from pyopencl.tools import (  # noqa
+        pytest_generate_tests_for_pyopencl
+        as pytest_generate_tests)
+from meshmode.mesh.generation import (  # noqa
+        generate_icosahedron, generate_box_mesh)
+from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry
+from meshmode.mesh.refinement import Refiner
+
+import logging
+logger = logging.getLogger(__name__)
+
+from functools import partial
+
+
+def gen_blob_mesh():
+    from meshmode.mesh.io import generate_gmsh, FileSource
+    return generate_gmsh(
+            FileSource("blob-2d.step"), 2, order=1,
+            force_ambient_dim=2,
+            other_options=[
+                "-string", "Mesh.CharacteristicLengthMax = %s;" % 0.2]
+            )
+
+
+def random_refine_flags(fract, mesh):
+    all_els = list(range(mesh.nelements))
+
+    flags = np.zeros(mesh.nelements)
+    from random import shuffle
+    shuffle(all_els)
+    for i in range(int(mesh.nelements * fract)):
+        flags[all_els[i]] = 1
+
+    return flags
+
+
+def uniform_refine_flags(mesh):
+    return np.ones(mesh.nelements)
+
+
+@pytest.mark.parametrize(("case_name", "mesh_gen", "flag_gen", "num_generations"), [
+    # Fails?
+    # ("icosahedron",
+    #     partial(generate_icosahedron, 1, order=1),
+    #     partial(random_refine_flags, 0.4),
+    #     3),
+
+    ("rect2d_rand",
+        partial(generate_box_mesh, (
+            np.linspace(0, 1, 3),
+            np.linspace(0, 1, 3),
+            ), order=1),
+        partial(random_refine_flags, 0.4),
+        4),
+
+    ("rect2d_unif",
+        partial(generate_box_mesh, (
+            np.linspace(0, 1, 2),
+            np.linspace(0, 1, 2),
+            ), order=1),
+        uniform_refine_flags,
+        3),
+
+    ("blob2d_rand",
+        gen_blob_mesh,
+        partial(random_refine_flags, 0.4),
+        4),
+
+    ("rect3d_rand",
+        partial(generate_box_mesh, (
+            np.linspace(0, 1, 2),
+            np.linspace(0, 1, 3),
+            np.linspace(0, 1, 2),
+            ), order=1),
+        partial(random_refine_flags, 0.4),
+        3),
+
+    ("rect3d_unif",
+        partial(generate_box_mesh, (
+            np.linspace(0, 1, 2),
+            np.linspace(0, 1, 2)), order=1),
+        uniform_refine_flags,
+        3),
+    ])
+def test_refinement(case_name, mesh_gen, flag_gen, num_generations):
+    from random import seed
+    seed(13)
+
+    mesh = mesh_gen()
+
+    r = Refiner(mesh)
+
+    for igen in range(num_generations):
+        flags = flag_gen(mesh)
+        mesh = r.refine(flags)
+
+        check_nodal_adj_against_geometry(mesh)
+
+
+if __name__ == "__main__":
+    import sys
+    if len(sys.argv) > 1:
+        exec(sys.argv[1])
+    else:
+        from py.test.cmdline import main
+        main([__file__])
+
+# vim: fdm=marker