diff --git a/meshmode/mesh/tools.py b/meshmode/mesh/tools.py
new file mode 100644
index 0000000000000000000000000000000000000000..a67a0fb8a310c0b7626a66159b567ae85762759c
--- /dev/null
+++ b/meshmode/mesh/tools.py
@@ -0,0 +1,48 @@
+from __future__ import division
+from __future__ import absolute_import
+from six.moves import range
+
+__copyright__ = "Copyright (C) 2010,2012,2013 Andreas Kloeckner, Michael Tom"
+
+__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 numpy as np
+from pytools.spatial_btree import SpatialBinaryTreeBucket
+
+
+def make_element_lookup_tree(mesh, eps=1e-12):
+    from meshmode.mesh.processing import find_bounding_box
+    bbox_min, bbox_max = find_bounding_box(mesh)
+    bbox_min -= eps
+    bbox_max += eps
+
+    tree = SpatialBinaryTreeBucket(bbox_min, bbox_max)
+
+    for igrp, grp in enumerate(mesh.groups):
+        for iel_grp in range(grp.nelements):
+            el_vertices = mesh.vertices[:, grp.vertex_indices[iel_grp]]
+
+            el_bbox_min = np.min(el_vertices, axis=-1) - eps
+            el_bbox_max = np.max(el_vertices, axis=-1) + eps
+
+            tree.insert((igrp, iel_grp), (el_bbox_min, el_bbox_max))
+
+    return tree
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index e7d3e94115a300fe087fb3d847324a9d3e74d9c5..9455098e9d1acf1f037463317db2ceb03bea35be 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -377,6 +377,29 @@ def test_as_python():
     assert mesh == mesh_2
 
 
+def test_lookup_tree(do_plot=False):
+    from meshmode.mesh.generation import make_curve_mesh, cloverleaf
+    mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, 1000), order=3)
+
+    from meshmode.mesh.tools import make_element_lookup_tree
+    tree = make_element_lookup_tree(mesh)
+
+    from meshmode.mesh.processing import find_bounding_box
+    bbox_min, bbox_max = find_bounding_box(mesh)
+
+    extent = bbox_max-bbox_min
+
+    for i in range(20):
+        pt = bbox_min + np.random.rand(2) * extent
+        print(pt)
+        for igrp, iel in tree.generate_matches(pt):
+            print(igrp, iel)
+
+    if do_plot:
+        with open("tree.dat", "w") as outf:
+            tree.visualize(outf)
+
+
 if __name__ == "__main__":
     import sys
     if len(sys.argv) > 1: