diff --git a/doc/conf.py b/doc/conf.py
index 866653625acb09a79319be6bec7c07872081f38f..5d9f499e4c3a0cc10b43ec14611fe88365bde549 100644
--- a/doc/conf.py
+++ b/doc/conf.py
@@ -280,5 +280,6 @@ intersphinx_mapping = {
     'https://documen.tician.de/pyopencl': None,
     'https://documen.tician.de/meshpy': None,
     'https://documen.tician.de/modepy': None,
-    'https://documen.tician.de/loopy': None
+    'https://documen.tician.de/loopy': None,
+    'https://tisaac.gitlab.io/recursivenodes/': None,
 }
diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py
index d0f3823c0b9bb368fa215a7187f39594bb190def..e317e31130237ad2adf5173af4a7c0d23d0f2b9d 100644
--- a/meshmode/discretization/poly_element.py
+++ b/meshmode/discretization/poly_element.py
@@ -40,6 +40,7 @@ Group types
 .. autoclass:: InterpolatoryQuadratureSimplexElementGroup
 .. autoclass:: QuadratureSimplexElementGroup
 .. autoclass:: PolynomialWarpAndBlendElementGroup
+.. autoclass:: PolynomialRecursiveNodesElementGroup
 .. autoclass:: PolynomialEquidistantSimplexElementGroup
 .. autoclass:: LegendreGaussLobattoTensorProductElementGroup
 
@@ -49,6 +50,7 @@ Group factories
 .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory
 .. autoclass:: QuadratureSimplexGroupFactory
 .. autoclass:: PolynomialWarpAndBlendGroupFactory
+.. autoclass:: PolynomialRecursiveNodesGroupFactory
 .. autoclass:: PolynomialEquidistantSimplexGroupFactory
 .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory
 """
@@ -207,6 +209,8 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup):
     polynomials in :math:`P^k`, hence usable for differentiation and
     interpolation. Interpolation nodes edge-clustered for avoidance of Runge
     phenomena. Nodes are present on the boundary of the simplex.
+
+    Uses :func:`modepy.warp_and_blend_nodes`.
     """
     @property
     @memoize_method
@@ -223,6 +227,41 @@ class PolynomialWarpAndBlendElementGroup(_MassMatrixQuadratureElementGroup):
         return result
 
 
+class PolynomialRecursiveNodesElementGroup(_MassMatrixQuadratureElementGroup):
+    """Elemental discretization with a number of nodes matching the number of
+    polynomials in :math:`P^k`, hence usable for differentiation and
+    interpolation. Interpolation nodes edge-clustered for avoidance of Runge
+    phenomena. Nodes are present on the boundary of the simplex.
+
+    Supports a choice of the base *family* of 1D nodes, see the documentation
+    of the *family* argument to :func:`recursivenodes.recursive_nodes`.
+
+    Requires :mod:`recursivenodes` to be installed.
+
+    .. [Isaac20] Tobin Isaac. Recursive, parameter-free, explicitly defined
+        interpolation nodes for simplices.
+        `Arxiv preprint <https://arxiv.org/abs/2002.09421>`__.
+
+    .. versionadded:: 2020.2
+    """
+    def __init__(self, mesh_el_group, order, family, index):
+        super().__init__(mesh_el_group, order, index)
+        self.family = family
+
+    @property
+    @memoize_method
+    def unit_nodes(self):
+        dim = self.mesh_el_group.dim
+
+        from recursivenodes import recursive_nodes
+        result = recursive_nodes(dim, self.order, self.family,
+                domain="biunit").T.copy()
+
+        dim2, nunit_nodes = result.shape
+        assert dim2 == dim
+        return result
+
+
 class PolynomialEquidistantSimplexElementGroup(_MassMatrixQuadratureElementGroup):
     """Elemental discretization with a number of nodes matching the number of
     polynomials in :math:`P^k`, hence usable for differentiation and
@@ -340,6 +379,20 @@ class PolynomialWarpAndBlendGroupFactory(HomogeneousOrderBasedGroupFactory):
     group_class = PolynomialWarpAndBlendElementGroup
 
 
+class PolynomialRecursiveNodesGroupFactory(HomogeneousOrderBasedGroupFactory):
+    def __init__(self, order, family):
+        self.order = order
+        self.family = family
+
+    def __call__(self, mesh_el_group, index):
+        if not isinstance(mesh_el_group, _MeshSimplexElementGroup):
+            raise TypeError("only mesh element groups of type '%s' "
+                    "are supported" % _MeshSimplexElementGroup.__name__)
+
+        return PolynomialRecursiveNodesElementGroup(
+                mesh_el_group, self.order, self.family, index)
+
+
 class PolynomialEquidistantSimplexGroupFactory(HomogeneousOrderBasedGroupFactory):
     """
     .. versionadded:: 2016.1
diff --git a/requirements.txt b/requirements.txt
index 4c3ff3392cea23471a0778e46ae7625bd6b7094e..574e50949645b761b6f3579c87fbb4e1dbfecc5e 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,4 +1,5 @@
 numpy
+recursivenodes
 git+https://gitlab.tiker.net/inducer/pytools.git#egg=pytools
 git+https://gitlab.tiker.net/inducer/gmsh_interop.git#egg=gmsh_interop
 git+https://gitlab.tiker.net/inducer/modepy.git#egg=modepy
diff --git a/setup.py b/setup.py
index 810d2560b674251f4b6b72d617e4fba74f38c692..1de30e57b763e1c01e91b5660e71023f93845238 100644
--- a/setup.py
+++ b/setup.py
@@ -44,6 +44,7 @@ def main():
               "numpy",
               "modepy",
               "gmsh_interop",
+              "recursivenodes",
               "six",
               "pytools>=2020.3",
               "pytest>=2.3",
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 17773c2086e0419705fac109175d77ebb34c0220..3417d0aa5128f7df3eb22f5c912a675dfabcb79e 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -22,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+from functools import partial
 from six.moves import range
 import numpy as np
 import numpy.linalg as la
@@ -35,6 +36,7 @@ from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup
 from meshmode.discretization.poly_element import (
         InterpolatoryQuadratureSimplexGroupFactory,
         PolynomialWarpAndBlendGroupFactory,
+        PolynomialRecursiveNodesGroupFactory,
         PolynomialEquidistantSimplexGroupFactory,
         )
 from meshmode.mesh import Mesh, BTAG_ALL
@@ -207,7 +209,9 @@ def test_box_boundary_tags(dim, nelem, group_factory):
 
 @pytest.mark.parametrize("group_factory", [
     InterpolatoryQuadratureSimplexGroupFactory,
-    PolynomialWarpAndBlendGroupFactory
+    PolynomialWarpAndBlendGroupFactory,
+    partial(PolynomialRecursiveNodesGroupFactory, family="lgl"),
+    #partial(PolynomialRecursiveNodesGroupFactory, family="gc"),
     ])
 @pytest.mark.parametrize("boundary_tag", [
     BTAG_ALL,
@@ -300,7 +304,7 @@ def test_boundary_interpolation(ctx_factory, group_factory, boundary_tag,
     print(eoc_rec)
     assert (
             eoc_rec.order_estimate() >= order-order_slack
-            or eoc_rec.max_error() < 1e-14)
+            or eoc_rec.max_error() < 3e-14)
 
 # }}}