diff --git a/examples/plot-connectivity.py b/examples/plot-connectivity.py
index 0faaa9b3319aca1345da4b507b4eb31ab6ef78b4..a96fc6f704b02c0c37398e05ca17f9c54cc0e1d5 100644
--- a/examples/plot-connectivity.py
+++ b/examples/plot-connectivity.py
@@ -18,11 +18,11 @@ def main():
     mesh = generate_icosahedron(1, order=order)
     #mesh = generate_torus(3, 1, order=order)
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.poly_element import \
             PolynomialWarpAndBlendGroupFactory
 
-    discr = Discretization(
+    discr = NodalDiscretization(
             actx, mesh, PolynomialWarpAndBlendGroupFactory(order))
 
     from meshmode.discretization.visualization import make_visualizer
diff --git a/examples/simple-dg.py b/examples/simple-dg.py
index e1d63ef69d19cbb8741c23fd9934cb1670bd4f60..cc04ef3f095168a55b4f2ff2ced566c0ef8a9d19 100644
--- a/examples/simple-dg.py
+++ b/examples/simple-dg.py
@@ -59,11 +59,11 @@ class DGDiscretization:
     def __init__(self, actx, mesh, order):
         self.order = order
 
-        from meshmode.discretization import Discretization
+        from meshmode.discretization.nodal import NodalDiscretization
         from meshmode.discretization.poly_element import \
                 PolynomialWarpAndBlendGroupFactory
         self.group_factory = PolynomialWarpAndBlendGroupFactory(order=order)
-        self.volume_discr = Discretization(actx, mesh, self.group_factory)
+        self.volume_discr = NodalDiscretization(actx, mesh, self.group_factory)
 
         assert self.volume_discr.dim == 2
 
diff --git a/examples/to_firedrake.py b/examples/to_firedrake.py
index 395211ca77f736e70d795a928d17988f3e4ca32c..0a961c9840395317cb227842502d198b06e5f2e5 100644
--- a/examples/to_firedrake.py
+++ b/examples/to_firedrake.py
@@ -61,11 +61,11 @@ def main():
 
     order = 3
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.poly_element import \
         InterpolatoryQuadratureSimplexGroupFactory
     group_factory = InterpolatoryQuadratureSimplexGroupFactory(order=order)
-    discr = Discretization(actx, mesh, group_factory)
+    discr = NodalDiscretization(actx, mesh, group_factory)
 
     # Get our solution: we will use
     # Real(e^z) = Real(e^{x+iy})
diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index dc949ebb85f096c5483c8e30a1a1af2cb0703369..6b9b92dc052e49ab2e617ae9c23a352bbfcbb773 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -1,4 +1,7 @@
-__copyright__ = "Copyright (C) 2013-2020 Andreas Kloeckner"
+__copyright__ = """
+Copyright (C) 2013-2021 Andreas Kloeckner
+Copyright (C) 2021 Thomas Gibson
+"""
 
 __license__ = """
 Permission is hereby granted, free of charge, to any person obtaining a copy
@@ -23,7 +26,6 @@ THE SOFTWARE.
 import numpy as np
 
 from pytools import memoize_in, memoize_method
-from pytools.obj_array import make_obj_array
 from meshmode.array_context import ArrayContext, make_loopy_program
 
 # underscored because it shouldn't be imported from here.
@@ -32,7 +34,7 @@ from meshmode.dof_array import DOFArray as _DOFArray
 __doc__ = """
 .. autoclass:: ElementGroupBase
 .. autoclass:: InterpolatoryElementGroupBase
-.. autoclass:: Discretization
+.. autoclass:: DiscretizationBase
 """
 
 
@@ -161,8 +163,8 @@ class InterpolatoryElementGroupBase(ElementGroupBase):
 # }}}
 
 
-class Discretization:
-    """An unstructured composite discretization.
+class DiscretizationBase:
+    """An unstructured composite discretization base class.
 
     .. attribute:: real_dtype
 
@@ -184,8 +186,6 @@ class Discretization:
     .. automethod:: empty_like
     .. automethod:: zeros_like
 
-    .. automethod:: nodes()
-
     .. automethod:: num_reference_derivative
 
     .. automethod:: quad_weights
@@ -344,50 +344,9 @@ class Discretization:
                         )["result"])
                 for grp in self.groups))
 
-    @memoize_method
-    def nodes(self):
-        r"""
-        :returns: object array of shape ``(ambient_dim,)`` containing
-            :class:`~meshmode.dof_array.DOFArray`\ s of node coordinates.
-        """
-
-        actx = self._setup_actx
-
-        @memoize_in(actx, (Discretization, "nodes_prg"))
-        def prg():
-            return make_loopy_program(
-                """{[iel,idof,j]:
-                    0<=iel<nelements and
-                    0<=idof<ndiscr_nodes and
-                    0<=j<nmesh_nodes}""",
-                """
-                    result[iel, idof] = \
-                        sum(j, resampling_mat[idof, j] * nodes[iel, j])
-                    """,
-                name="nodes")
-
-        def resample_mesh_nodes(grp, iaxis):
-            # TODO: would be nice to have the mesh use an array context already
-            nodes = actx.from_numpy(grp.mesh_el_group.nodes[iaxis])
-
-            grp_unit_nodes = grp.unit_nodes.reshape(-1)
-            meg_unit_nodes = grp.mesh_el_group.unit_nodes.reshape(-1)
-
-            tol = 10 * np.finfo(grp_unit_nodes.dtype).eps
-            if (grp_unit_nodes.shape == meg_unit_nodes.shape
-                    and np.linalg.norm(grp_unit_nodes - meg_unit_nodes) < tol):
-                return nodes
-
-            return actx.call_loopy(
-                    prg(),
-                    resampling_mat=actx.from_numpy(grp.from_mesh_interp_matrix()),
-                    nodes=nodes,
-                    )["result"]
 
-        return make_obj_array([
-            _DOFArray(None, tuple([
-                actx.freeze(resample_mesh_nodes(grp, iaxis)) for grp in self.groups
-                ]))
-            for iaxis in range(self.ambient_dim)])
+# For backwards compatibility, we need to be sure we export the
+# alias for NodalDiscretization
+from meshmode.discretization.nodal import Discretization
 
 # vim: fdm=marker
diff --git a/meshmode/discretization/modal.py b/meshmode/discretization/modal.py
new file mode 100644
index 0000000000000000000000000000000000000000..00dc7a34ece9dce3edbbb4497bb288578b1279bd
--- /dev/null
+++ b/meshmode/discretization/modal.py
@@ -0,0 +1,72 @@
+__copyright__ = """
+Copyright (C) 2013-2021 Andreas Kloeckner
+Copyright (C) 2021 Thomas Gibson
+"""
+
+__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.
+"""
+
+
+from pytools import memoize_in, memoize_method
+from pytools.obj_array import make_obj_array
+from meshmode.array_context import make_loopy_program
+
+# underscored because it shouldn't be imported from here.
+from meshmode.dof_array import DOFArray as _DOFArray
+
+from meshmode.discretization import DiscretizationBase
+
+
+class ModalDiscretization(DiscretizationBase):
+    """An unstructured composite modal discretization."""
+
+    @memoize_method
+    def modes(self):
+        r"""
+        :returns: object array of shape ``(ambient_dim,)`` containing
+            :class:`~meshmode.dof_array.DOFArray`\ s of modal coefficients.
+        """
+
+        actx = self._setup_actx
+
+        @memoize_in(actx, (ModalDiscretization, "modes_prg"))
+        def prg():
+            return make_loopy_program(
+                """{[iel,idof,j]:
+                    0<=iel<nelements and
+                    0<=idof<ndiscr_nodes and
+                    0<=j<nmesh_nodes}""",
+                """
+                    result[iel, idof] = \
+                        sum(j, resampling_mat[idof, j] * modes[iel, j])
+                    """,
+                name="modes")
+
+        return make_obj_array([
+            _DOFArray(None, tuple(
+                actx.freeze(
+                    actx.call_loopy(
+                        prg(),
+                        resampling_mat=actx.from_numpy(
+                            grp.from_mesh_interp_matrix()),
+                        nodes=actx.from_numpy(grp.mesh_el_group.nodes[iaxis])
+                        )["result"])
+                for grp in self.groups))
+            for iaxis in range(self.ambient_dim)])
diff --git a/meshmode/discretization/nodal.py b/meshmode/discretization/nodal.py
new file mode 100644
index 0000000000000000000000000000000000000000..366e15a1c1add26d15445f17bdb7956db9c8f0f5
--- /dev/null
+++ b/meshmode/discretization/nodal.py
@@ -0,0 +1,83 @@
+__copyright__ = """
+Copyright (C) 2013-2021 Andreas Kloeckner
+Copyright (C) 2021 Thomas Gibson
+"""
+
+__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.
+"""
+
+
+from pytools import memoize_in, memoize_method
+from pytools.obj_array import make_obj_array
+from meshmode.array_context import make_loopy_program
+
+# underscored because it shouldn't be imported from here.
+from meshmode.dof_array import DOFArray as _DOFArray
+
+from meshmode.discretization import DiscretizationBase
+
+
+class NodalDiscretization(DiscretizationBase):
+    """An unstructured composite nodal discretization.
+
+    .. automethod:: nodes()
+    """
+
+    @memoize_method
+    def nodes(self):
+        r"""
+        :returns: object array of shape ``(ambient_dim,)`` containing
+            :class:`~meshmode.dof_array.DOFArray`\ s of node coordinates.
+        """
+
+        actx = self._setup_actx
+
+        @memoize_in(actx, (NodalDiscretization, "nodes_prg"))
+        def prg():
+            return make_loopy_program(
+                """{[iel,idof,j]:
+                    0<=iel<nelements and
+                    0<=idof<ndiscr_nodes and
+                    0<=j<nmesh_nodes}""",
+                """
+                    result[iel, idof] = \
+                        sum(j, resampling_mat[idof, j] * nodes[iel, j])
+                    """,
+                name="nodes")
+
+        return make_obj_array([
+            _DOFArray(None, tuple(
+                actx.freeze(
+                    actx.call_loopy(
+                        prg(),
+                        resampling_mat=actx.from_numpy(
+                            grp.from_mesh_interp_matrix()),
+                        nodes=actx.from_numpy(grp.mesh_el_group.nodes[iaxis])
+                        )["result"])
+                for grp in self.groups))
+            for iaxis in range(self.ambient_dim)])
+
+
+# Create a temporary alias for NodalDiscretization until completely
+# phased out in libraries using meshmode.
+def Discretization(*args, **kwargs):  # noqa: N802
+    from warnings import warn
+    warn("Discretization will be depreciated in favor of NodalDiscretization")
+    return NodalDiscretization(*args, **kwargs)
diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py
index 0e438d9da16af0800819e7ac5fffc27fccbbd946..6d739b155453ba85777f312411f768be7ebaa2f9 100644
--- a/meshmode/interop/nodal_dg.py
+++ b/meshmode/interop/nodal_dg.py
@@ -102,7 +102,7 @@ class NodalDGContext:
 
         self.octave.eval(f"StartUp{mesh.dim}D;")
 
-    def get_discr(self, actx) -> meshmode.discretization.Discretization:
+    def get_discr(self, actx) -> meshmode.discretization.nodal.NodalDiscretization:
         """Get a discretization with nodes exactly matching the ones used
         by the nodal-DG code.
 
@@ -144,10 +144,10 @@ class NodalDGContext:
 
         mesh = Mesh(vertices=vertices, groups=[egroup], is_conforming=True)
 
-        from meshmode.discretization import Discretization
+        from meshmode.discretization.nodal import NodalDiscretization
         from meshmode.discretization.poly_element import (
                 PolynomialGivenNodesGroupFactory)
-        return Discretization(actx, mesh,
+        return NodalDiscretization(actx, mesh,
                 PolynomialGivenNodesGroupFactory(order, unit_nodes))
 
     def push_dof_array(self, name, ary: meshmode.dof_array.DOFArray):
diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py
index a60486268cc4f5a811fed8b888fab66a5fe2af7b..20bb70b5d01b21ecf4beed5f47ae4e60552d5ad7 100644
--- a/meshmode/mesh/visualization.py
+++ b/meshmode/mesh/visualization.py
@@ -309,8 +309,8 @@ def vtk_visualize_mesh(actx, mesh, filename, vtk_high_order=True):
 
     from meshmode.discretization.poly_element import \
             PolynomialWarpAndBlendGroupFactory
-    from meshmode.discretization import Discretization
-    discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+    from meshmode.discretization.nodal import NodalDiscretization as NDiscr
+    discr = NDiscr(actx, mesh, PolynomialWarpAndBlendGroupFactory(order))
 
     from meshmode.discretization.visualization import make_visualizer
     vis = make_visualizer(actx, discr, order, force_equidistant=vtk_high_order)
diff --git a/test/test_array.py b/test/test_array.py
index 2f089a05b6ff99813707730c6a23b37743ed8f8f..b6c157d1e447aa0593c0e7bc12b83137111c5ec9 100644
--- a/test/test_array.py
+++ b/test/test_array.py
@@ -27,7 +27,7 @@ from meshmode.array_context import (  # noqa
         pytest_generate_tests_for_pyopencl_array_context
         as pytest_generate_tests)
 
-from meshmode.discretization import Discretization
+from meshmode.discretization.nodal import NodalDiscretization
 from meshmode.discretization.poly_element import PolynomialWarpAndBlendGroupFactory
 from meshmode.dof_array import flatten, unflatten, DOFArray
 
@@ -44,7 +44,7 @@ def test_array_context_np_workalike(actx_factory):
     mesh = generate_regular_rect_mesh(
             a=(-0.5,)*2, b=(0.5,)*2, n=(8,)*2, order=3)
 
-    discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3))
+    discr = NodalDiscretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3))
 
     for sym_name, n_args in [
             ("sin", 1),
@@ -88,7 +88,7 @@ def test_dof_array_arithmetic_same_as_numpy(actx_factory):
     mesh = generate_regular_rect_mesh(
             a=(-0.5,)*2, b=(0.5,)*2, n=(3,)*2, order=1)
 
-    discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3))
+    discr = NodalDiscretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(3))
 
     def get_real(ary):
         return ary.real
diff --git a/test/test_chained.py b/test/test_chained.py
index 70296c6f561666f2b777ed4c0d82a6518b593e40..dea34243aa0437ab6a672ae88c6332120df6436f 100644
--- a/test/test_chained.py
+++ b/test/test_chained.py
@@ -70,10 +70,10 @@ def create_discretization(actx, ndim,
         raise ValueError(f"unsupported dimension: {ndim}")
 
     # create discretization
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.poly_element import \
             InterpolatoryQuadratureSimplexGroupFactory
-    discr = Discretization(actx, mesh,
+    discr = NodalDiscretization(actx, mesh,
             InterpolatoryQuadratureSimplexGroupFactory(order))
 
     return discr
diff --git a/test/test_firedrake_interop.py b/test/test_firedrake_interop.py
index 045978c89f115f04d702750bd8f02d493699b1bb..c4313783a5b7f091f911d090e10b0d9b922ade1e 100644
--- a/test/test_firedrake_interop.py
+++ b/test/test_firedrake_interop.py
@@ -29,7 +29,7 @@ from pyopencl.tools import (  # noqa
 
 from meshmode.array_context import PyOpenCLArrayContext
 
-from meshmode.discretization import Discretization
+from meshmode.discretization.nodal import NodalDiscretization
 from meshmode.discretization.poly_element import (
     InterpolatoryQuadratureSimplexGroupFactory)
 
@@ -192,7 +192,7 @@ def test_to_fd_consistency(ctx_factory, mm_mesh, fspace_degree):
     actx = PyOpenCLArrayContext(queue)
 
     factory = InterpolatoryQuadratureSimplexGroupFactory(fspace_degree)
-    discr = Discretization(actx, mm_mesh, factory)
+    discr = NodalDiscretization(actx, mm_mesh, factory)
     fdrake_connection = build_connection_to_firedrake(discr)
     fdrake_fspace = fdrake_connection.firedrake_fspace()
     # Check consistency
@@ -538,7 +538,7 @@ def test_to_fd_transfer(ctx_factory, fspace_degree, mesh_name, mesh_pars, dim):
 
         # Make discr and connect it to firedrake
         factory = InterpolatoryQuadratureSimplexGroupFactory(fspace_degree)
-        discr = Discretization(actx, mm_mesh, factory)
+        discr = NodalDiscretization(actx, mm_mesh, factory)
 
         fdrake_connection = build_connection_to_firedrake(discr)
         fdrake_fspace = fdrake_connection.firedrake_fspace()
@@ -658,7 +658,7 @@ def test_to_fd_idempotency(ctx_factory, mm_mesh, fspace_degree):
 
     # Make a function space and a function with unique values at each node
     factory = InterpolatoryQuadratureSimplexGroupFactory(fspace_degree)
-    discr = Discretization(actx, mm_mesh, factory)
+    discr = NodalDiscretization(actx, mm_mesh, factory)
     fdrake_connection = build_connection_to_firedrake(discr)
     fdrake_mesh = fdrake_connection.firedrake_fspace().mesh()
     dtype = fdrake_mesh.coordinates.dat.data.dtype
diff --git a/test/test_mesh.py b/test/test_mesh.py
index 1bf7685188afcaf4baf257b4a1b6065a96958061..28c9b02610478ddc57d686870d861b7f3fa9307c 100644
--- a/test/test_mesh.py
+++ b/test/test_mesh.py
@@ -42,10 +42,10 @@ def test_nonequal_rect_mesh_generation(actx_factory, visualize=False):
     mesh = mgen.generate_regular_rect_mesh(
             a=(0, 0)*2, b=(5, 3), n=(10, 6,), order=3)
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.poly_element import \
             PolynomialWarpAndBlendGroupFactory as GroupFactory
-    discr = Discretization(actx, mesh, GroupFactory(3))
+    discr = NodalDiscretization(actx, mesh, GroupFactory(3))
 
     if visualize:
         from meshmode.discretization.visualization import make_visualizer
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 1a50edfc71ef0413140ae2bb7263f2646fb1e2da..bc3b73f8c2e9aa6e1102cce84c9c92585b9e0901 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -112,8 +112,8 @@ def test_parallel_vtk_file(actx_factory, dim):
     else:
         raise ValueError("unknown dimensionality")
 
-    from meshmode.discretization import Discretization
-    discr = Discretization(actx, mesh,
+    from meshmode.discretization.nodal import NodalDiscretization
+    discr = NodalDiscretization(actx, mesh,
             InterpolatoryQuadratureSimplexGroupFactory(target_order))
 
     from meshmode.discretization.visualization import make_visualizer
@@ -193,8 +193,8 @@ def test_visualizers(actx_factory, dim, group_cls):
     else:
         group_factory = LegendreGaussLobattoTensorProductGroupFactory
 
-    from meshmode.discretization import Discretization
-    discr = Discretization(actx, mesh, group_factory(target_order))
+    from meshmode.discretization.nodal import NodalDiscretization
+    discr = NodalDiscretization(actx, mesh, group_factory(target_order))
 
     nodes = thaw(actx, discr.nodes())
     f = actx.np.sqrt(sum(nodes**2)) + 1j*nodes[0]
@@ -419,7 +419,7 @@ def test_boundary_interpolation(actx_factory, group_factory, boundary_tag,
     else:
         group_cls = SimplexElementGroup
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.connection import (
             make_face_restriction, check_connection)
 
@@ -468,7 +468,7 @@ def test_boundary_interpolation(actx_factory, group_factory, boundary_tag,
 
         # }}}
 
-        vol_discr = Discretization(actx, mesh, group_factory(order))
+        vol_discr = NodalDiscretization(actx, mesh, group_factory(order))
         print("h=%s -> %d elements" % (
                 h, sum(mgrp.nelements for mgrp in mesh.groups)))
 
@@ -533,7 +533,7 @@ def test_all_faces_interpolation(actx_factory, group_factory,
     else:
         group_cls = SimplexElementGroup
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.connection import (
             make_face_restriction, make_face_to_all_faces_embedding,
             check_connection)
@@ -574,7 +574,7 @@ def test_all_faces_interpolation(actx_factory, group_factory,
 
         # }}}
 
-        vol_discr = Discretization(actx, mesh, group_factory(order))
+        vol_discr = NodalDiscretization(actx, mesh, group_factory(order))
         print("h=%s -> %d elements" % (
                 h, sum(mgrp.nelements for mgrp in mesh.groups)))
 
@@ -656,7 +656,7 @@ def test_opposite_face_interpolation(actx_factory, group_factory,
     else:
         group_cls = SimplexElementGroup
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.connection import (
             make_face_restriction, make_opposite_face_connection,
             check_connection)
@@ -705,7 +705,7 @@ def test_opposite_face_interpolation(actx_factory, group_factory,
 
         # }}}
 
-        vol_discr = Discretization(actx, mesh, group_factory(order))
+        vol_discr = NodalDiscretization(actx, mesh, group_factory(order))
         print("h=%s -> %d elements" % (
                 h, sum(mgrp.nelements for mgrp in mesh.groups)))
 
@@ -783,8 +783,8 @@ def test_orientation_3d(actx_factory, what, mesh_gen_func, visualize=False):
 
     logger.info("%d elements", mesh.nelements)
 
-    from meshmode.discretization import Discretization
-    discr = Discretization(actx, mesh,
+    from meshmode.discretization.nodal import NodalDiscretization
+    discr = NodalDiscretization(actx, mesh,
             PolynomialWarpAndBlendGroupFactory(3))
 
     from pytential import bind, sym
@@ -865,9 +865,9 @@ def test_merge_and_map(actx_factory, group_cls, visualize=False):
     mesh4 = mesh3.copy()
 
     if visualize:
-        from meshmode.discretization import Discretization
+        from meshmode.discretization.nodal import NodalDiscretization
         actx = actx_factory()
-        discr = Discretization(actx, mesh4, discr_grp_factory)
+        discr = NodalDiscretization(actx, mesh4, discr_grp_factory)
 
         from meshmode.discretization.visualization import make_visualizer
         vis = make_visualizer(actx, discr, 3, element_shrink_factor=0.8)
@@ -910,8 +910,8 @@ def test_sanity_single_element(actx_factory, dim, mesh_order, group_cls,
     mg = group_cls(mesh_order, vertex_indices, nodes, dim=dim)
     mesh = Mesh(vertices, [mg], is_conforming=True)
 
-    from meshmode.discretization import Discretization
-    vol_discr = Discretization(actx, mesh, group_factory)
+    from meshmode.discretization.nodal import NodalDiscretization
+    vol_discr = NodalDiscretization(actx, mesh, group_factory)
 
     # {{{ volume calculation check
 
@@ -985,10 +985,10 @@ def test_sanity_qhull_nd(actx_factory, dim, order):
     mesh = from_vertices_and_simplices(dtri.points.T, dtri.simplices,
             fix_orientation=True)
 
-    from meshmode.discretization import Discretization
-    low_discr = Discretization(actx, mesh,
+    from meshmode.discretization.nodal import NodalDiscretization
+    low_discr = NodalDiscretization(actx, mesh,
             PolynomialEquidistantSimplexGroupFactory(order))
-    high_discr = Discretization(actx, mesh,
+    high_discr = NodalDiscretization(actx, mesh,
             PolynomialEquidistantSimplexGroupFactory(order+1))
 
     from meshmode.discretization.connection import make_same_mesh_connection
@@ -1050,8 +1050,8 @@ def test_sanity_balls(actx_factory, src_file, dim, mesh_order, visualize=False):
 
         # {{{ discretizations and connections
 
-        from meshmode.discretization import Discretization
-        vol_discr = Discretization(actx, mesh,
+        from meshmode.discretization.nodal import NodalDiscretization
+        vol_discr = NodalDiscretization(actx, mesh,
                 InterpolatoryQuadratureSimplexGroupFactory(quad_order))
 
         from meshmode.discretization.connection import make_face_restriction
@@ -1151,10 +1151,12 @@ def test_box_mesh(actx_factory, visualize=False):
     mesh = mgen.generate_box_mesh(3*(np.linspace(0, 1, 5),))
 
     if visualize:
-        from meshmode.discretization import Discretization
+        from meshmode.discretization.nodal import NodalDiscretization
+        from meshmode.discretization.poly_element import \
+                PolynomialWarpAndBlendGroupFactory
 
         actx = actx_factory()
-        discr = Discretization(actx, mesh,
+        discr = NodalDiscretization(actx, mesh,
                 PolynomialWarpAndBlendGroupFactory(7))
 
         from meshmode.discretization.visualization import make_visualizer
@@ -1404,8 +1406,8 @@ def test_vtk_overwrite(actx_factory):
 
     mesh = mgen.generate_torus(10.0, 2.0, order=target_order)
 
-    from meshmode.discretization import Discretization
-    discr = Discretization(
+    from meshmode.discretization.nodal import NodalDiscretization
+    discr = NodalDiscretization(
             actx, mesh,
             InterpolatoryQuadratureSimplexGroupFactory(target_order))
 
@@ -1478,8 +1480,8 @@ def test_mesh_without_vertices(actx_factory):
     mesh = refine_uniformly(mesh, 1)
 
     # make sure the world doesn't end
-    from meshmode.discretization import Discretization
-    discr = Discretization(actx, mesh,
+    from meshmode.discretization.nodal import NodalDiscretization
+    discr = NodalDiscretization(actx, mesh,
             InterpolatoryQuadratureSimplexGroupFactory(4))
     thaw(actx, discr.nodes())
 
@@ -1605,8 +1607,8 @@ def test_mesh_multiple_groups(actx_factory, ambient_dim, visualize=False):
         import matplotlib.pyplot as plt
         plt.savefig("test_mesh_multiple_groups_2d_elements.png", dpi=300)
 
-    from meshmode.discretization import Discretization
-    discr = Discretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order))
+    from meshmode.discretization.nodal import NodalDiscretization
+    discr = NodalDiscretization(actx, mesh, PolynomialWarpAndBlendGroupFactory(order))
 
     if visualize:
         group_id = discr.empty(actx, dtype=np.int32)
@@ -1696,9 +1698,9 @@ def test_mesh_with_interior_unit_nodes(actx_factory, ambient_dim):
     assert mesh.facial_adjacency_groups
     assert mesh.nodal_adjacency
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.poly_element import QuadratureSimplexGroupFactory
-    discr = Discretization(actx, mesh,
+    discr = NodalDiscretization(actx, mesh,
             QuadratureSimplexGroupFactory(order))
 
     from meshmode.discretization.connection import make_face_restriction
diff --git a/test/test_partition.py b/test/test_partition.py
index 77db5abecb852c80cbd5d774bb46bd537f2f806f..4232349e1a8c3ece0a7856976b19222dc1ecaf11 100644
--- a/test/test_partition.py
+++ b/test/test_partition.py
@@ -103,8 +103,8 @@ def test_partition_interpolation(actx_factory, dim, mesh_pars,
             for i_remote_part in neighbors:
                 connected_parts.add((i_local_part, i_remote_part))
 
-        from meshmode.discretization import Discretization
-        vol_discrs = [Discretization(actx, part_meshes[i], group_factory(order))
+        from meshmode.discretization.nodal import NodalDiscretization
+        vol_discrs = [NodalDiscretization(actx, part_meshes[i], group_factory(order))
                         for i in range(num_parts)]
 
         from meshmode.mesh import BTAG_PARTITION
@@ -348,8 +348,8 @@ def _test_mpi_boundary_swap(dim, order, num_groups):
     queue = cl.CommandQueue(cl_ctx)
     actx = PyOpenCLArrayContext(queue)
 
-    from meshmode.discretization import Discretization
-    vol_discr = Discretization(actx, local_mesh, group_factory)
+    from meshmode.discretization.nodal import NodalDiscretization
+    vol_discr = NodalDiscretization(actx, local_mesh, group_factory)
 
     from meshmode.distributed import get_connected_partitions
     connected_parts = get_connected_partitions(local_mesh)
diff --git a/test/test_refinement.py b/test/test_refinement.py
index ca748e4f29ee911a8f4934131e2cec7602ad0f73..42ebf9b99382e9555cca1f9db8dc37df88fea106 100644
--- a/test/test_refinement.py
+++ b/test/test_refinement.py
@@ -202,7 +202,7 @@ def test_refinement_connection(
     # discretization order
     order = 5
 
-    from meshmode.discretization import Discretization
+    from meshmode.discretization.nodal import NodalDiscretization
     from meshmode.discretization.connection import (
             make_refinement_connection, check_connection)
 
@@ -249,7 +249,7 @@ def test_refinement_connection(
 
             return result
 
-        discr = Discretization(actx, mesh, group_factory(order))
+        discr = NodalDiscretization(actx, mesh, group_factory(order))
 
         refiner = refiner_cls(mesh)
         flags = refine_flags(mesh)