diff --git a/setup.py b/setup.py
index 19924d68d83aa0f8dbd087810e6ac3d1108b8628..c1e5f698e535fd1891cdedf6ff6e00c6765dd1e2 100644
--- a/setup.py
+++ b/setup.py
@@ -44,6 +44,7 @@ setup(name="sumpy",
           "pytools>=2013.5.6",
           "boxtree>=2013.1",
           "pytest>=2.3",
+          "pyfmmlib>=2016.1",
           "six",
 
           # If this causes issues, see:
diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py
new file mode 100644
index 0000000000000000000000000000000000000000..b0e4ef26538fddff804b0cf48266f80734f7d9e0
--- /dev/null
+++ b/sumpy/expansion/level_to_order.py
@@ -0,0 +1,93 @@
+from __future__ import division
+
+__copyright__ = "Copyright (C) 2016 Matt Wala"
+
+__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.
+"""
+
+__doc__ = """
+.. autofunction:: h2d_level_to_order_lookup
+.. autofunction:: l2d_level_to_order_lookup
+"""
+
+import numpy as np
+
+
+def h2d_level_to_order_lookup(tree, helmholtz_k, epsilon):
+    """
+    Compute a mapping from level number to expansion order,
+    Helmholtz 2D case.
+
+    This wraps the function *h2dterms* from :mod:`pyfmmlib`.
+
+    :arg tree: An instance of :class:`boxtree.Tree`.
+    :arg helmholtz_k: Helmholtz parameter
+    :arg epsilon: Precision
+
+    :return: A :class:`numpy.array` of length `tree.nlevels`
+    """
+
+    if tree.dimensions != 2:
+        raise ValueError("tree must be 2d")
+
+    orders = np.empty(tree.nlevels, dtype=int)
+    bbox_area = np.max(
+        tree.bounding_box[1] - tree.bounding_box[0]) ** 2
+
+    from pyfmmlib import h2dterms
+    for level in range(tree.nlevels):
+        nterms, ier = h2dterms(bbox_area / 2 ** level, helmholtz_k, epsilon)
+        if ier != 0:
+            raise RuntimeError(
+                "h2dterms returned error code {ier}".format(ier=ier))
+        orders[level] = nterms
+
+    return orders
+
+
+def l2d_level_to_order_lookup(tree, epsilon):
+    """
+    Compute a mapping from level number to expansion order,
+    Laplace 2D case.
+
+    This wraps the function *l2dterms* from :mod:`pyfmmlib`.
+
+    :arg tree: An instance of :class:`boxtree.Tree`.
+    :arg epsilon: Precision
+
+    :return: A :class:`numpy.array` of length `tree.nlevels`
+    """
+
+    if tree.dimensions != 2:
+        raise ValueError("tree must be 2d")
+
+    from pyfmmlib import l2dterms
+    nterms, ier = l2dterms(epsilon)
+    if ier != 0:
+        raise RuntimeError(
+            "l2dterms returned error code {ier}".format(ier=ier))
+
+    orders = np.empty(tree.nlevels, dtype=int)
+    orders.fill(nterms)
+
+    return orders
+
+
+# vim: fdm=marker
diff --git a/test/test_fmm.py b/test/test_fmm.py
index ce26622a3e0d6e2130e2f30944c82bf2f55a4962..fa99c56ebe08c8ad8566e7658b98a783fac5bc3b 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -36,6 +36,10 @@ from sumpy.expansion.multipole import (
 from sumpy.expansion.local import (
     VolumeTaylorLocalExpansion, H2DLocalExpansion)
 
+from sumpy.expansion.level_to_order import (
+    h2d_level_to_order_lookup,
+    l2d_level_to_order_lookup)
+
 import pytest
 
 import logging
@@ -50,6 +54,34 @@ else:
     faulthandler.enable()
 
 
+@pytest.mark.parametrize("lookup_func, extra_args", [
+    (h2d_level_to_order_lookup, (100,)),
+    (l2d_level_to_order_lookup, ()),
+    ])
+def test_level_to_order_lookup(ctx_getter, lookup_func, extra_args):
+    ctx = ctx_getter()
+    queue = cl.CommandQueue(ctx)
+
+    nsources = 500
+    dtype = np.float64
+    epsilon = 1e-9
+
+    from boxtree.tools import (
+            make_normal_particle_array as p_normal)
+
+    sources = p_normal(queue, nsources, 2, dtype, seed=15)
+    from boxtree import TreeBuilder
+    tb = TreeBuilder(ctx)
+
+    tree, _ = tb(queue, sources, max_particles_in_box=30, debug=True)
+
+    levels = lookup_func(tree, *(extra_args + (epsilon,)))
+
+    assert len(levels) == tree.nlevels
+    # Should be monotonically nonincreasing.
+    assert all(levels[i] >= levels[i+1] for i in range(tree.nlevels - 1))
+
+
 @pytest.mark.parametrize("knl, local_expn_class, mpole_expn_class", [
     (LaplaceKernel(2), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),
     (LaplaceKernel(3), VolumeTaylorLocalExpansion, VolumeTaylorMultipoleExpansion),