diff --git a/sumpy/expansion/level_to_order.py b/sumpy/expansion/level_to_order.py
index 8f012e55827f7676e2429f76502f5146c085387c..5e29b86138ec705ef57747929dead3403160b17a 100644
--- a/sumpy/expansion/level_to_order.py
+++ b/sumpy/expansion/level_to_order.py
@@ -23,73 +23,61 @@ THE SOFTWARE.
 """
 
 __doc__ = """
-.. autofunction:: h2d_level_to_order_lookup
-.. autofunction:: l2d_level_to_order_lookup
-
+.. autoclass:: FMMLibExpansionOrderFinder
 .. autoclass:: SimpleExpansionOrderFinder
 """
 
 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`
+class FMMLibExpansionOrderFinder(object):
+    r"""Return expansion orders that meet the tolerance for a given level
+    using routines wrapped from :mod:`pyfmmlib`.
     """
 
-    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 __init__(self, tol, extra_order=0):
+        """
+        :arg tol: error tolerance
+        :arg extra_order: order increase to accommodate, say, the taking of
+            derivatives of the FMM expansions.
+        """
+        self.tol = tol
+        self.extra_order = extra_order
 
+    def __call__(self, kernel, kernel_args, tree, level):
+        import pyfmmlib
 
-def l2d_level_to_order_lookup(tree, epsilon):
-    """
-    Compute a mapping from level number to expansion order,
-    Laplace 2D case.
+        from sumpy.kernel import LaplaceKernel, HelmholtzKernel
 
-    This wraps the function *l2dterms* from :mod:`pyfmmlib`.
+        assert isinstance(kernel, (LaplaceKernel, HelmholtzKernel))
+        assert tree.dimensions in (2, 3)
 
-    :arg tree: An instance of :class:`boxtree.Tree`.
-    :arg epsilon: Precision
+        if isinstance(kernel, LaplaceKernel):
+            if tree.dimensions == 2:
+                nterms, ier = pyfmmlib.l2dterms(self.tol)
+                if ier:
+                    raise RuntimeError("l2dterms returned error code '%d'" % ier)
 
-    :return: A :class:`numpy.array` of length `tree.nlevels`
-    """
+            elif tree.dimensions == 3:
+                nterms, ier = pyfmmlib.l3dterms(self.tol)
+                if ier:
+                    raise RuntimeError("l3dterms returned error code '%d'" % ier)
 
-    if tree.dimensions != 2:
-        raise ValueError("tree must be 2d")
+        elif isinstance(kernel, HelmholtzKernel):
+            helmholtz_k = dict(kernel_args)[kernel.helmholtz_k_name]
+            size = tree.root_extent / 2 ** level
 
-    from pyfmmlib import l2dterms
-    nterms, ier = l2dterms(epsilon)
-    if ier != 0:
-        raise RuntimeError(
-            "l2dterms returned error code {ier}".format(ier=ier))
+            if tree.dimensions == 2:
+                nterms, ier = pyfmmlib.h2dterms(size, helmholtz_k, self.tol)
+                if ier:
+                    raise RuntimeError("h2dterms returned error code '%d'" % ier)
 
-    orders = np.empty(tree.nlevels, dtype=int)
-    orders.fill(nterms)
+            elif tree.dimensions == 3:
+                nterms, ier = pyfmmlib.h3dterms(size, helmholtz_k, self.tol)
+                if ier:
+                    raise RuntimeError("h3dterms returned error code '%d'" % ier)
 
-    return orders
+        return nterms + self.extra_order
 
 
 class SimpleExpansionOrderFinder(object):
@@ -114,7 +102,7 @@ class SimpleExpansionOrderFinder(object):
             extra_order=1):
         """
         :arg extra_order: order increase to accommodate, say, the taking of
-            oderivatives f the FMM expansions.
+            derivatives of the FMM expansions.
         """
         self.tol = tol
 
diff --git a/test/test_fmm.py b/test/test_fmm.py
index 71e3f044432b4bd612337bdd55d56e421540bcc4..2bc6d69d6bcafedf338cc78aea0591d6d5ac3845 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -42,10 +42,6 @@ from sumpy.expansion.local import (
     LaplaceConformingVolumeTaylorLocalExpansion,
     HelmholtzConformingVolumeTaylorLocalExpansion)
 
-from sumpy.expansion.level_to_order import (
-    h2d_level_to_order_lookup,
-    l2d_level_to_order_lookup)
-
 import pytest
 
 import logging
@@ -60,35 +56,6 @@ 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):
-    pytest.importorskip("pyfmmlib")
-    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(2), LaplaceConformingVolumeTaylorLocalExpansion,
diff --git a/test/test_misc.py b/test/test_misc.py
index 33e8535a8cbce064263958962421d9b43e77ee6c..a25b2e800ab6f9942395821cf0cb4345e80b6ccf 100644
--- a/test/test_misc.py
+++ b/test/test_misc.py
@@ -124,7 +124,9 @@ class FakeTree:
         self.stick_out_factor = stick_out_factor
 
 
-@pytest.mark.parametrize("knl", [LaplaceKernel(2), HelmholtzKernel(2)])
+@pytest.mark.parametrize("knl", [
+        LaplaceKernel(2), HelmholtzKernel(2),
+        LaplaceKernel(3), HelmholtzKernel(3)])
 def test_order_finder(knl):
     from sumpy.expansion.level_to_order import SimpleExpansionOrderFinder
 
@@ -136,6 +138,28 @@ def test_order_finder(knl):
         for level in range(30)]
     print(orders)
 
+    # Order should not increase with level
+    assert (np.diff(orders) <= 0).all()
+
+
+@pytest.mark.parametrize("knl", [
+        LaplaceKernel(2), HelmholtzKernel(2),
+        LaplaceKernel(3), HelmholtzKernel(3)])
+def test_fmmlib_order_finder(knl):
+    pytest.importorskip("pyfmmlib")
+    from sumpy.expansion.level_to_order import FMMLibExpansionOrderFinder
+
+    ofind = FMMLibExpansionOrderFinder(1e-5)
+
+    tree = FakeTree(knl.dim, 200, 0.5)
+    orders = [
+        ofind(knl, frozenset([("k", 5)]), tree, level)
+        for level in range(30)]
+    print(orders)
+
+    # Order should not increase with level
+    assert (np.diff(orders) <= 0).all()
+
 
 # {{{ expansion toys p2e2e2p test cases