diff --git a/examples/simple-dg.py b/examples/simple-dg.py
index 324f2a925..ac109b6be 100644
--- a/examples/simple-dg.py
+++ b/examples/simple-dg.py
@@ -273,12 +273,13 @@ def get_local_face_mass_matrix(self, afgrp, volgrp, dtype):
                     afgrp.nunit_dofs),
                 dtype=dtype)
 
-        from modepy.tools import UNIT_VERTICES
         import modepy as mp
-        for iface, fvi in enumerate(
-                volgrp.mesh_el_group.face_vertex_indices()):
-            face_vertices = UNIT_VERTICES[volgrp.dim][np.array(fvi)].T
-            matrix[:, iface, :] = mp.nodal_face_mass_matrix(
+        shape = mp.Simplex(volgrp.dim)
+        unit_vertices = mp.unit_vertices_for_shape(shape).T
+
+        for face in mp.faces_for_shape(shape):
+            face_vertices = unit_vertices[np.array(face.volume_vertex_indices)].T
+            matrix[:, face.face_index, :] = mp.nodal_face_mass_matrix(
                     volgrp.basis(), volgrp.unit_nodes, afgrp.unit_nodes,
                     volgrp.order,
                     face_vertices)
diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py
index f164d5d72..7ec816b23 100644
--- a/meshmode/discretization/__init__.py
+++ b/meshmode/discretization/__init__.py
@@ -91,15 +91,15 @@ def nelements(self):
 
     @property
     def nunit_dofs(self):
-        """The number of (for now: nodal) degrees of freedom ("DOFs") associated with a
-        single element.
+        """The number of (for now: nodal) degrees of freedom ("DOFs")
+        associated with a single element.
         """
         return self.unit_nodes.shape[-1]
 
     @property
     def ndofs(self):
-        """The total number of (for now: nodal) degrees of freedom ("DOFs") associated with
-        the element group.
+        """The total number of (for now: nodal) degrees of freedom ("DOFs")
+        associated with the element group.
         """
         return self.nunit_dofs * self.nelements
 
diff --git a/meshmode/discretization/connection/face.py b/meshmode/discretization/connection/face.py
index 9eddc07c7..e3494c019 100644
--- a/meshmode/discretization/connection/face.py
+++ b/meshmode/discretization/connection/face.py
@@ -74,8 +74,7 @@ def _build_boundary_connection(actx, vol_discr, bdry_discr, connection_data,
             bdry_grp = bdry_discr.groups[ibdry_grp]
             data = connection_data[igrp, face_id]
 
-            bdry_unit_nodes_01 = (bdry_grp.unit_nodes + 1)*0.5
-            result_unit_nodes = (np.dot(data.A, bdry_unit_nodes_01).T + data.b).T
+            result_unit_nodes = data.face.map_to_volume(bdry_grp.unit_nodes)
 
             batches.append(
                 InterpolationBatch(
@@ -140,9 +139,9 @@ def _get_face_vertices(mesh, boundary_tag):
 
         # }}}
     else:
-        # For FRESTR_INTERIOR_FACES, this is likely every vertex in the book.
+        # For FACE_RESTR_INTERIOR, this is likely every vertex in the book.
         # Don't ever bother trying to cut the list down.
-        # For FRESTR_ALL_FACES, it literally is every single vertex.
+        # For FACE_RESTR_ALL, it literally is every single vertex.
 
         return np.arange(mesh.nvertices, dtype=np.intp)
 
@@ -208,7 +207,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag,
 
     # }}}
 
-    from meshmode.mesh import Mesh, SimplexElementGroup
+    from meshmode.mesh import Mesh, _ModepyElementGroup
     bdry_mesh_groups = []
     connection_data = {}
 
@@ -220,9 +219,10 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag,
 
         mgrp = grp.mesh_el_group
 
-        if not isinstance(mgrp, SimplexElementGroup):
+        if not isinstance(mgrp, _ModepyElementGroup):
             raise NotImplementedError("can only take boundary of "
-                    "SimplexElementGroup-based meshes")
+                    "meshes based on SimplexElementGroup and "
+                    "TensorProductElementGroup")
 
         # {{{ pull together per-group face lists
 
@@ -257,44 +257,42 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag,
 
         # }}}
 
-        grp_face_vertex_indices = mgrp.face_vertex_indices()
-        grp_vertex_unit_coordinates = mgrp.vertex_unit_coordinates()
-
         batch_base = 0
 
-        # group by face_id
+        # group by face_index
 
-        for face_id in range(mgrp.nfaces):
-            batch_boundary_el_numbers_in_grp = np.array(
-                    [
-                        ibface_el
-                        for ibface_el, ibface_face in group_boundary_faces
-                        if ibface_face == face_id],
-                    dtype=np.intp)
+        for face in mgrp._modepy_faces:
+            batch_boundary_el_numbers_in_grp = np.array([
+                ibface_el
+                for ibface_el, ibface_face in group_boundary_faces
+                if ibface_face == face.face_index
+                ], dtype=np.intp)
 
             # {{{ preallocate arrays for mesh group
 
             nbatch_elements = len(batch_boundary_el_numbers_in_grp)
 
-            if per_face_groups or face_id == 0:
+            if per_face_groups or face.face_index == 0:
                 if per_face_groups:
                     ngroup_bdry_elements = nbatch_elements
                 else:
                     ngroup_bdry_elements = len(group_boundary_faces)
 
+                # make up some not-terrible nodes for the boundary Mesh
+                space = mp.space_for_shape(face, mgrp.order)
+                bdry_unit_nodes = mp.edge_clustered_nodes_for_space(space, face)
+
+                vol_basis = mp.basis_for_space(
+                        mgrp._modepy_space, mgrp._modepy_shape).functions
+
                 vertex_indices = np.empty(
-                        (ngroup_bdry_elements, mgrp.dim+1-1),
+                        (ngroup_bdry_elements, face.nvertices),
                         mgrp.vertex_indices.dtype)
 
-                bdry_unit_nodes = mp.warp_and_blend_nodes(mgrp.dim-1, mgrp.order)
-                bdry_unit_nodes_01 = (bdry_unit_nodes + 1)*0.5
-
-                vol_basis = mp.simplex_onb(mgrp.dim, mgrp.order)
-                nbdry_unit_nodes = bdry_unit_nodes_01.shape[-1]
+                nbdry_unit_nodes = bdry_unit_nodes.shape[-1]
                 nodes = np.empty(
                         (discr.ambient_dim, ngroup_bdry_elements, nbdry_unit_nodes),
                         dtype=np.float64)
-
             # }}}
 
             new_el_numbers = batch_base + np.arange(nbatch_elements)
@@ -303,24 +301,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag,
 
             # {{{ no per-element axes in these computations
 
-            # Find boundary vertex indices
-            loc_face_vertices = list(grp_face_vertex_indices[face_id])
-
-            # Find unit nodes for boundary element
-            face_vertex_unit_coordinates = \
-                    grp_vertex_unit_coordinates[loc_face_vertices]
-
-            # Find A, b such that A [e_1 e_2] + b = [r_1 r_2]
-            # (Notation assumes that the volume is 3D and the face is 2D.
-            # Code does not.)
-
-            b = face_vertex_unit_coordinates[0]
-            A = (  # noqa
-                    face_vertex_unit_coordinates[1:]
-                    - face_vertex_unit_coordinates[0]).T
-
-            face_unit_nodes = (np.dot(A, bdry_unit_nodes_01).T + b).T
-
+            face_unit_nodes = face.map_to_volume(bdry_unit_nodes)
             resampling_mat = mp.resampling_matrix(
                     vol_basis,
                     face_unit_nodes, mgrp.unit_nodes)
@@ -331,7 +312,7 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag,
 
             # Find vertex_indices
             glob_face_vertices = mgrp.vertex_indices[
-                    batch_boundary_el_numbers_in_grp][:, loc_face_vertices]
+                    batch_boundary_el_numbers_in_grp][:, face.volume_vertex_indices]
             vertex_indices[new_el_numbers] = \
                     vol_to_bdry_vertices[glob_face_vertices]
 
@@ -343,17 +324,16 @@ def make_face_restriction(actx, discr, group_factory, boundary_tag,
 
             # }}}
 
-            connection_data[igrp, face_id] = _ConnectionBatchData(
+            connection_data[igrp, face.face_index] = _ConnectionBatchData(
                     group_source_element_indices=batch_boundary_el_numbers_in_grp,
                     group_target_element_indices=new_el_numbers,
-                    A=A,
-                    b=b,
+                    face=face,
                     )
 
-            is_last_face = face_id + 1 == mgrp.nfaces
+            is_last_face = face.face_index + 1 == mgrp.nfaces
 
             if per_face_groups or is_last_face:
-                bdry_mesh_group = SimplexElementGroup(
+                bdry_mesh_group = type(mgrp)(
                         mgrp.order, vertex_indices, nodes,
                         unit_nodes=bdry_unit_nodes)
                 bdry_mesh_groups.append(bdry_mesh_group)
diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py
index 22cd60cd4..5714ca8df 100644
--- a/meshmode/discretization/poly_element.py
+++ b/meshmode/discretization/poly_element.py
@@ -41,6 +41,8 @@
 .. autoclass:: PolynomialRecursiveNodesElementGroup
 .. autoclass:: PolynomialEquidistantSimplexElementGroup
 .. autoclass:: PolynomialGivenNodesElementGroup
+
+.. autoclass:: GaussLegendreTensorProductElementGroup
 .. autoclass:: LegendreGaussLobattoTensorProductElementGroup
 .. autoclass:: EquidistantTensorProductElementGroup
 
@@ -48,12 +50,16 @@
 ^^^^^^^^^^^^^^^
 
 .. autoclass:: ElementGroupFactory
+.. autoclass:: OrderAndTypeBasedGroupFactory
+
 .. autoclass:: InterpolatoryQuadratureSimplexGroupFactory
 .. autoclass:: QuadratureSimplexGroupFactory
 .. autoclass:: PolynomialWarpAndBlendGroupFactory
 .. autoclass:: PolynomialRecursiveNodesGroupFactory
 .. autoclass:: PolynomialEquidistantSimplexGroupFactory
 .. autoclass:: PolynomialGivenNodesGroupFactory
+
+.. autoclass:: GaussLegendreTensorProductGroupFactory
 .. autoclass:: LegendreGaussLobattoTensorProductGroupFactory
 """
 
@@ -98,11 +104,22 @@ def diff_matrices(self):
 # {{{ concrete element groups for simplices
 
 class SimplexElementGroupBase(ElementGroupBase):
+    @property
+    @memoize_method
+    def _shape(self):
+        return mp.Simplex(self.dim)
+
+    @property
+    @memoize_method
+    def _space(self):
+        return mp.PN(self.dim, self.order)
+
     @memoize_method
     def from_mesh_interp_matrix(self):
         meg = self.mesh_el_group
+        meg_space = mp.PN(meg.dim, meg.order)
         return mp.resampling_matrix(
-                mp.simplex_best_available_basis(meg.dim, meg.order),
+                mp.basis_for_space(meg_space, self._shape).functions,
                 self.unit_nodes,
                 meg.unit_nodes)
 
@@ -112,28 +129,19 @@ class PolynomialSimplexElementGroupBase(PolynomialElementGroupBase,
     def is_orthogonal_basis(self):
         return self.dim <= 3
 
+    @property
     @memoize_method
-    def _mode_ids_and_basis(self):
-        # for now, see https://gitlab.tiker.net/inducer/modepy/-/merge_requests/14
-        import modepy.modes as modes
-        if self.dim <= 3:
-            return modes.simplex_onb_with_mode_ids(self.dim, self.order)
-        else:
-            return modes.simplex_monomial_basis_with_mode_ids(self.dim, self.order)
+    def _basis(self):
+        return mp.basis_for_space(self._space, self._shape)
 
     def basis(self):
-        mode_ids, basis = self._mode_ids_and_basis()
-        return basis
+        return self._basis.functions
 
     def mode_ids(self):
-        mode_ids, basis = self._mode_ids_and_basis()
-        return mode_ids
+        return self._basis.mode_ids
 
     def grad_basis(self):
-        if self.dim <= 3:
-            return mp.grad_simplex_onb(self.dim, self.order)
-        else:
-            return mp.grad_simplex_monomial_basis(self.dim, self.order)
+        return self._basis.gradients
 
 
 class InterpolatoryQuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase):
@@ -340,52 +348,153 @@ def unit_nodes(self):
 
         return self._unit_nodes
 
-
 # }}}
 
 
 # {{{ concrete element groups for tensor product elements
 
-class LegendreGaussLobattoTensorProductElementGroup(PolynomialElementGroupBase):
-    def is_orthogonal_basis(self):
-        return True
-
-    def basis(self):
-        from modepy.modes import tensor_product_basis, jacobi
-        from functools import partial
-        return tensor_product_basis(
-                self.dim, tuple(
-                    partial(jacobi, 0, 0, i)
-                    for i in range(self.order+1)))
-
-    def grad_basis(self):
-        raise NotImplementedError()
+class HypercubeElementGroupBase(ElementGroupBase):
+    @property
+    @memoize_method
+    def _shape(self):
+        return mp.Hypercube(self.dim)
 
     @property
     @memoize_method
-    def unit_nodes(self):
-        from modepy.nodes import tensor_product_nodes
-        from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
-        return tensor_product_nodes(
-                self.dim, legendre_gauss_lobatto_nodes(self.order))
+    def _space(self):
+        return mp.QN(self.dim, self.order)
 
     @memoize_method
     def from_mesh_interp_matrix(self):
         meg = self.mesh_el_group
+        meg_space = mp.QN(meg.dim, meg.order)
         return mp.resampling_matrix(
-                self.basis(),
+                mp.basis_for_space(meg_space, self._shape).functions,
                 self.unit_nodes,
                 meg.unit_nodes)
 
 
-class EquidistantTensorProductElementGroup(
-        LegendreGaussLobattoTensorProductElementGroup):
-    @property
+class TensorProductElementGroupBase(PolynomialElementGroupBase,
+        HypercubeElementGroupBase):
+    def __init__(self, mesh_el_group, order, index, *, basis, unit_nodes):
+        """
+        :arg basis: a :class:`modepy.TensorProductBasis`.
+        :arg unit_nodes: unit nodes for the tensor product, obtained by
+            using :func:`modepy.tensor_product_nodes`, for example.
+        """
+        super().__init__(mesh_el_group, order, index)
+
+        if basis._dim != mesh_el_group.dim:
+            raise ValueError("basis dimension does not match element group: "
+                    f"expected {mesh_el_group.dim}, got {basis._dim}.")
+
+        if unit_nodes.shape[0] != mesh_el_group.dim:
+            raise ValueError("unit node dimension does not match element group: "
+                    f"expected {mesh_el_group.dim}, got {unit_nodes.shape[0]}.")
+
+        self._basis = basis
+        self._unit_nodes = unit_nodes
+
+    def is_orthogonal_basis(self):
+        try:
+            # NOTE: meshmode kind of assumes that the basis is orthonormal
+            # with weight 1, which is why this check is stricter than expected.
+            return self._basis.orthonormality_weight() == 1
+        except mp.BasisNotOrthonormal:
+            return False
+
     @memoize_method
+    def mode_ids(self):
+        return self._basis.mode_ids
+
+    @memoize_method
+    def basis(self):
+        return self._basis.functions
+
+    @memoize_method
+    def grad_basis(self):
+        return self._basis.gradients
+
+    @property
     def unit_nodes(self):
-        from modepy.nodes import tensor_product_nodes, equidistant_nodes
-        return tensor_product_nodes(
-                self.dim, equidistant_nodes(1, self.order)[0])
+        return self._unit_nodes
+
+    @property
+    @memoize_method
+    def weights(self):
+        return np.dot(
+                self.mass_matrix(),
+                np.ones(len(self.basis())))
+
+
+class LegendreTensorProductElementGroup(TensorProductElementGroupBase):
+    def __init__(self, mesh_el_group, order, index, *, unit_nodes):
+        basis = mp.orthonormal_basis_for_space(
+                mp.QN(mesh_el_group.dim, order),
+                mp.Hypercube(mesh_el_group.dim))
+
+        super().__init__(mesh_el_group, order, index,
+                basis=basis,
+                unit_nodes=unit_nodes)
+
+
+class GaussLegendreTensorProductElementGroup(LegendreTensorProductElementGroup):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes matching the number of polynomials in the tensor
+    product basis, hence usable for differentiation and interpolation.
+
+    No interpolation nodes are present on the boundary of the hypercube.
+    """
+
+    def __init__(self, mesh_el_group, order, index):
+        self._quadrature_rule = mp.LegendreGaussTensorProductQuadrature(
+                order, mesh_el_group.dim)
+
+        super().__init__(mesh_el_group, order, index,
+                unit_nodes=self._quadrature_rule.nodes)
+
+    @property
+    def weights(self):
+        return self._quadrature_rule.weights
+
+
+class LegendreGaussLobattoTensorProductElementGroup(
+        LegendreTensorProductElementGroup):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes matching the number of polynomials in the tensor
+    product basis, hence usable for differentiation and interpolation.
+    Nodes sufficient for unisolvency are present on the boundary of the hypercube.
+
+    Uses :func:`~modepy.quadrature.jacobi_gauss.legendre_gauss_lobatto_nodes`.
+    """
+
+    def __init__(self, mesh_el_group, order, index):
+        from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
+        unit_nodes_1d = legendre_gauss_lobatto_nodes(order)
+
+        super().__init__(mesh_el_group, order, index,
+                unit_nodes=mp.tensor_product_nodes(
+                    [unit_nodes_1d] * mesh_el_group.dim)
+                )
+
+
+class EquidistantTensorProductElementGroup(LegendreTensorProductElementGroup):
+    """Elemental discretization supplying a high-order quadrature rule
+    with a number of nodes matching the number of polynomials in the tensor
+    product basis, hence usable for differentiation and interpolation.
+    Nodes sufficient for unisolvency are present on the boundary of the hypercube.
+
+    Uses :func:`~modepy.equidistant_nodes`.
+    """
+
+    def __init__(self, mesh_el_group, order, index):
+        from modepy.nodes import equidistant_nodes
+        unit_nodes_1d = equidistant_nodes(1, order)[0]
+
+        super().__init__(mesh_el_group, order, index,
+                unit_nodes=mp.tensor_product_nodes(
+                    [unit_nodes_1d] * mesh_el_group.dim)
+                )
 
 # }}}
 
@@ -483,9 +592,17 @@ def __call__(self, mesh_el_group, index):
 
         return PolynomialGivenNodesElementGroup(
                 mesh_el_group, self.order, self.unit_nodes, index)
+
 # }}}
 
 
+# {{{ group factories for tensor products
+
+class GaussLegendreTensorProductGroupFactory(HomogeneousOrderBasedGroupFactory):
+    mesh_group_class = _MeshTensorProductElementGroup
+    group_class = GaussLegendreTensorProductElementGroup
+
+
 class LegendreGaussLobattoTensorProductGroupFactory(
         HomogeneousOrderBasedGroupFactory):
     mesh_group_class = _MeshTensorProductElementGroup
@@ -493,5 +610,7 @@ class LegendreGaussLobattoTensorProductGroupFactory(
 
 # }}}
 
+# }}}
+
 
 # vim: fdm=marker
diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py
index 2426b25cd..1b2643ace 100644
--- a/meshmode/discretization/visualization.py
+++ b/meshmode/discretization/visualization.py
@@ -20,10 +20,14 @@
 THE SOFTWARE.
 """
 
+from functools import singledispatch
+
 import numpy as np
+
 from pytools import memoize_method, Record
 from meshmode.dof_array import flatten, thaw
 
+from modepy.shapes import Shape, Simplex, Hypercube
 
 __doc__ = """
 
@@ -117,6 +121,56 @@ def primitive_element_size(self):
 # }}}
 
 
+# {{{ vtk submeshes
+
+@singledispatch
+def vtk_submesh_for_shape(shape: Shape, node_tuples):
+    raise NotImplementedError(type(shape).__name__)
+
+
+@vtk_submesh_for_shape.register(Simplex)
+def _(shape: Simplex, node_tuples):
+    import modepy as mp
+    return mp.submesh_for_shape(shape, node_tuples)
+
+
+@vtk_submesh_for_shape.register(Hypercube)
+def _(shape: Hypercube, node_tuples):
+    node_tuple_to_index = {nt: i for i, nt in enumerate(node_tuples)}
+
+    # NOTE: this can't use mp.submesh_for_shape because VTK vertex order is
+    # counterclockwise instead of z order
+    el_offsets = {
+            1: [(0,), (1,)],
+            2: [(0, 0), (1, 0), (1, 1), (0, 1)],
+            3: [
+                (0, 0, 0),
+                (1, 0, 0),
+                (1, 1, 0),
+                (0, 1, 0),
+                (0, 0, 1),
+                (1, 0, 1),
+                (1, 1, 1),
+                (0, 1, 1),
+                ]
+            }[shape.dim]
+
+    from pytools import add_tuples
+    elements = []
+    for origin in node_tuples:
+        try:
+            elements.append(tuple(
+                node_tuple_to_index[add_tuples(origin, offset)]
+                for offset in el_offsets
+                ))
+        except KeyError:
+            pass
+
+    return elements
+
+# }}}
+
+
 # {{{ vtk connectivity
 
 class VTKConnectivity:
@@ -155,51 +209,24 @@ def tensor_cell_types(self):
                 }
 
     def connectivity_for_element_group(self, grp):
-        from pytools import (
-                generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam,
-                generate_nonnegative_integer_tuples_below as gnitb)
-        from meshmode.mesh import TensorProductElementGroup, SimplexElementGroup
+        import modepy as mp
+        from meshmode.mesh import _ModepyElementGroup
 
-        if isinstance(grp.mesh_el_group, SimplexElementGroup):
-            node_tuples = list(gnitstam(grp.order, grp.dim))
+        if isinstance(grp.mesh_el_group, _ModepyElementGroup):
+            shape = grp.mesh_el_group._modepy_shape
+            space = type(grp.mesh_el_group._modepy_space)(grp.dim, grp.order)
+            node_tuples = mp.node_tuples_for_space(space)
 
-            from modepy.tools import simplex_submesh
             el_connectivity = np.array(
-                    simplex_submesh(node_tuples),
+                    vtk_submesh_for_shape(shape, node_tuples),
                     dtype=np.intp)
 
-            vtk_cell_type = self.simplex_cell_types[grp.dim]
-
-        elif isinstance(grp.mesh_el_group, TensorProductElementGroup):
-            node_tuples = list(gnitb(grp.order+1, grp.dim))
-            node_tuple_to_index = {
-                    nt: i for i, nt in enumerate(node_tuples)}
-
-            def add_tuple(a, b):
-                return tuple(ai+bi for ai, bi in zip(a, b))
-
-            el_offsets = {
-                    1: [(0,), (1,)],
-                    2: [(0, 0), (1, 0), (1, 1), (0, 1)],
-                    3: [
-                        (0, 0, 0),
-                        (1, 0, 0),
-                        (1, 1, 0),
-                        (0, 1, 0),
-                        (0, 0, 1),
-                        (1, 0, 1),
-                        (1, 1, 1),
-                        (0, 1, 1),
-                        ]
-                    }[grp.dim]
-
-            el_connectivity = np.array([
-                    [
-                        node_tuple_to_index[add_tuple(origin, offset)]
-                        for offset in el_offsets]
-                    for origin in gnitb(grp.order, grp.dim)])
-
-            vtk_cell_type = self.tensor_cell_types[grp.dim]
+            if isinstance(shape, Simplex):
+                vtk_cell_type = self.simplex_cell_types[shape.dim]
+            elif isinstance(shape, Hypercube):
+                vtk_cell_type = self.tensor_cell_types[shape.dim]
+            else:
+                raise TypeError(f"unsupported shape: {type(shape)}")
 
         else:
             raise NotImplementedError("visualization for element groups "
diff --git a/meshmode/interop/nodal_dg.py b/meshmode/interop/nodal_dg.py
index 1bfca3f6f..5f43732cc 100644
--- a/meshmode/interop/nodal_dg.py
+++ b/meshmode/interop/nodal_dg.py
@@ -113,11 +113,11 @@ def get_discr(self, actx) -> meshmode.discretization.Discretization:
             unit_nodes_arrays = self.octave.eval(
                     f"Nodes{dim}D(N)", nout=dim, verbose=False)
 
-            equilat_to_biunit_func_name = (
+            equilat_to_unit_func_name = (
                     "".join(self.AXES[:dim] + ["to"] + self.REF_AXES[:dim]))
 
             unit_nodes_arrays = self.octave.feval(
-                    equilat_to_biunit_func_name, *unit_nodes_arrays,
+                    equilat_to_unit_func_name, *unit_nodes_arrays,
                     nout=dim, verbose=False)
 
             unit_nodes = np.array([a.reshape(-1) for a in unit_nodes_arrays])
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index ddb4f4160..ea298b185 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -199,32 +199,33 @@ class MeshElementGroup(Record):
 
     def __init__(self, order, vertex_indices, nodes,
             element_nr_base=None, node_nr_base=None,
-            unit_nodes=None, dim=None):
+            unit_nodes=None, dim=None, **kwargs):
         """
         :arg order: the maximum total degree used for interpolation.
-        :arg nodes: ``[ambient_dim, nelements, nunit_nodes]``
+        :arg nodes: ``(ambient_dim, nelements, nunit_nodes)``
             The nodes are assumed to be mapped versions of *unit_nodes*.
-        :arg unit_nodes: ``[dim, nunit_nodes]``
-            The unit nodes of which *nodes* is a mapped
-            version.
+        :arg unit_nodes: ``(dim, nunit_nodes)`` The unit nodes of which *nodes*
+            is a mapped version.
 
         Do not supply *element_nr_base* and *node_nr_base*, they will be
         automatically assigned.
         """
 
-        Record.__init__(self,
+        super().__init__(
             order=order,
             vertex_indices=vertex_indices,
             nodes=nodes,
             unit_nodes=unit_nodes,
-            element_nr_base=element_nr_base, node_nr_base=node_nr_base)
+            element_nr_base=element_nr_base, node_nr_base=node_nr_base,
+            **kwargs)
 
-    def copy(self, **kwargs):
+    def get_copy_kwargs(self, **kwargs):
         if "element_nr_base" not in kwargs:
             kwargs["element_nr_base"] = None
         if "node_nr_base" not in kwargs:
             kwargs["node_nr_base"] = None
-        return Record.copy(self, **kwargs)
+
+        return super().get_copy_kwargs(**kwargs)
 
     @property
     def dim(self):
@@ -255,20 +256,20 @@ def nunit_nodes(self):
 
     @property
     def is_affine(self):
-        raise NotImplementedError()
+        raise NotImplementedError
 
     def face_vertex_indices(self):
         """Return a tuple of tuples indicating which vertices
         (in mathematically positive ordering) make up each face
         of an element in this group.
         """
-        raise NotImplementedError()
+        raise NotImplementedError
 
     def vertex_unit_coordinates(self):
         """Return an array of shape ``(nfaces, dim)`` with the unit
         coordinates of each vertex.
         """
-        raise NotImplementedError()
+        raise NotImplementedError
 
     @property
     def nfaces(self):
@@ -290,134 +291,94 @@ def __ne__(self, other):
 # }}}
 
 
-# {{{ simplex
-
-class SimplexElementGroup(MeshElementGroup):
-    """
-    .. automethod:: __init__
-    """
+# {{{ modepy-based element group
 
+class _ModepyElementGroup(MeshElementGroup):
     def __init__(self, order, vertex_indices, nodes,
             element_nr_base=None, node_nr_base=None,
-            unit_nodes=None, dim=None):
+            unit_nodes=None, dim=None, **kwargs):
         """
         :arg order: the maximum total degree used for interpolation.
-        :arg nodes: ``[ambient_dim, nelements, nunit_nodes]``
+        :arg nodes: ``(ambient_dim, nelements, nunit_nodes)``
             The nodes are assumed to be mapped versions of *unit_nodes*.
-        :arg unit_nodes: ``[dim, nunit_nodes]``
-            The unit nodes of which *nodes* is a mapped
-            version. If unspecified, the nodes from
-            :func:`modepy.warp_and_blend_nodes` for *dim*
-            are assumed. These must be in unit coordinates
-            as defined in :mod:`modepy`.
-        :arg dim: only used if *unit_nodes* is None, to get
+        :arg unit_nodes: ``(dim, nunit_nodes)`` The unit nodes of which
+            *nodes* is a mapped version. If unspecified, the nodes from
+            :func:`modepy.edge_clustered_nodes_for_space` are assumed.
+            These must be in unit coordinates as defined in :mod:`modepy`.
+        :arg dim: only used if *unit_nodes* is *None*, to get
             the default unit nodes.
 
         Do not supply *element_nr_base* and *node_nr_base*, they will be
         automatically assigned.
         """
 
-        if unit_nodes is None:
+        if unit_nodes is not None:
+            _dim = unit_nodes.shape[0]
+            if dim is not None and _dim != dim:
+                raise ValueError("'dim' does not match 'unit_nodes' dimension")
+            else:
+                dim = _dim
+        else:
             if dim is None:
-                raise TypeError("'dim' must be passed "
-                        "if 'unit_nodes' is not passed")
+                raise TypeError("'dim' must be passed if 'unit_nodes' is not passed")
 
-            if dim <= 3:
-                unit_nodes = mp.warp_and_blend_nodes(dim, order)
-            else:
-                unit_nodes = mp.equidistant_nodes(dim, order)
+        # dim is now usable
+        shape = self._modepy_shape_cls(dim)
+        space = mp.space_for_shape(shape, order)
 
-        dims = unit_nodes.shape[0]
+        if unit_nodes is None:
+            unit_nodes = mp.edge_clustered_nodes_for_space(space, shape)
+
+        if nodes is not None:
+            if unit_nodes.shape[-1] != nodes.shape[-1]:
+                raise ValueError(
+                        "'nodes' has wrong number of unit nodes per element."
+                        f" expected {unit_nodes.shape[-1]}, "
+                        f" but got {nodes.shape[-1]}.")
 
         if vertex_indices is not None:
             if not issubclass(vertex_indices.dtype.type, np.integer):
-                raise TypeError("vertex_indices must be integral")
+                raise TypeError("'vertex_indices' must be integral")
 
-            if vertex_indices.shape[-1] != dims+1:
-                raise ValueError("vertex_indices has wrong number of vertices per "
-                        "element. expected: %d, got: %d" % (dims+1,
-                            vertex_indices.shape[-1]))
+            if vertex_indices.shape[-1] != shape.nvertices:
+                raise ValueError(
+                        "'vertex_indices' has wrong number of vertices per element."
+                        f" expected {shape.nvertices},"
+                        f" got {vertex_indices.shape[-1]}")
 
         super().__init__(order, vertex_indices, nodes,
-                element_nr_base, node_nr_base, unit_nodes, dim)
+                element_nr_base=element_nr_base,
+                node_nr_base=node_nr_base,
+                unit_nodes=unit_nodes,
+                dim=dim,
+                _modepy_shape=shape,
+                _modepy_space=space)
 
     @property
     @memoize_method
-    def is_affine(self):
-        return is_affine_simplex_group(self)
+    def _modepy_faces(self):
+        return mp.faces_for_shape(self._modepy_shape)
 
     def face_vertex_indices(self):
-        if self.dim == 1:
-            return (
-                (0,),
-                (1,),
-                )
-        elif self.dim == 2:
-            return (
-                (0, 1),
-                (2, 0),
-                (1, 2),
-                )
-        elif self.dim == 3:
-            return (
-                (0, 2, 1),
-                (0, 1, 3),
-                (0, 3, 2),
-                (1, 2, 3)
-                )
-        else:
-            raise NotImplementedError("dim=%d" % self.dim)
+        return tuple(face.volume_vertex_indices for face in self._modepy_faces)
 
     def vertex_unit_coordinates(self):
-        from modepy.tools import unit_vertices
-        return unit_vertices(self.dim)
+        return mp.unit_vertices_for_shape(self._modepy_shape).T
 
 # }}}
 
 
-# {{{ tensor-product
+class SimplexElementGroup(_ModepyElementGroup):
+    _modepy_shape_cls = mp.Simplex
 
-class TensorProductElementGroup(MeshElementGroup):
-    """
-    .. automethod:: __init__
-    """
-
-    def __init__(self, order, vertex_indices, nodes,
-            element_nr_base=None, node_nr_base=None,
-            unit_nodes=None):
-        """
-        :arg order: the maximum total degree used for interpolation.
-        :arg nodes: ``[ambient_dim, nelements, nunit_nodes]``
-            The nodes are assumed to be mapped versions of *unit_nodes*.
-        :arg unit_nodes: ``[dim, nunit_nodes]``
-            The unit nodes of which *nodes* is a mapped
-            version.
-
-        Do not supply *element_nr_base* and *node_nr_base*, they will be
-        automatically assigned.
-        """
-
-        dims = unit_nodes.shape[0]
-
-        if vertex_indices is not None:
-            if not issubclass(vertex_indices.dtype.type, np.integer):
-                raise TypeError("vertex_indices must be integral")
-
-            if vertex_indices.shape[-1] != 2**dims:
-                raise ValueError("vertex_indices has wrong number of vertices per "
-                        "element. expected: %d, got: %d" % (2**dims,
-                            vertex_indices.shape[-1]))
-
-        super().__init__(order, vertex_indices, nodes,
-                element_nr_base, node_nr_base, unit_nodes)
-
-    def face_vertex_indices(self):
-        raise NotImplementedError()
+    @property
+    @memoize_method
+    def is_affine(self):
+        return is_affine_simplex_group(self)
 
-    def vertex_unit_coordinates(self):
-        raise NotImplementedError()
 
-# }}}
+class TensorProductElementGroup(_ModepyElementGroup):
+    _modepy_shape_cls = mp.Hypercube
 
 # }}}
 
@@ -973,15 +934,20 @@ def __ne__(self, other):
 
 # {{{ node-vertex consistency test
 
-def _test_node_vertex_consistency_simplex(mesh, mgrp, tol):
+def _test_node_vertex_consistency_resampling(mesh, mgrp, tol):
     if mesh.vertices is None:
         return True
 
     if mgrp.nelements == 0:
         return True
 
+    if isinstance(mgrp, _ModepyElementGroup):
+        basis = mp.basis_for_space(mgrp._modepy_space, mgrp._modepy_shape).functions
+    else:
+        raise TypeError(f"unsupported group type: {type(mgrp).__name__}")
+
     resampling_mat = mp.resampling_matrix(
-            mp.simplex_best_available_basis(mgrp.dim, mgrp.order),
+            basis,
             mgrp.vertex_unit_coordinates().T,
             mgrp.unit_nodes)
 
@@ -1015,8 +981,8 @@ def _test_node_vertex_consistency(mesh, tol):
     """
 
     for mgrp in mesh.groups:
-        if isinstance(mgrp, SimplexElementGroup):
-            assert _test_node_vertex_consistency_simplex(mesh, mgrp, tol)
+        if isinstance(mgrp, _ModepyElementGroup):
+            assert _test_node_vertex_consistency_resampling(mesh, mgrp, tol)
         else:
             from warnings import warn
             warn("not implemented: node-vertex consistency check for '%s'"
@@ -1437,11 +1403,10 @@ def is_affine_simplex_group(group, abs_tol=None):
                 type(group).__name__)
 
     # get matrices
-    basis = mp.simplex_best_available_basis(group.dim, group.order)
-    grad_basis = mp.grad_simplex_best_available_basis(group.dim, group.order)
-
-    vinv = la.inv(mp.vandermonde(basis, group.unit_nodes))
-    diff = mp.differentiation_matrices(basis, grad_basis, group.unit_nodes)
+    basis = mp.basis_for_space(group._modepy_space, group._modepy_shape)
+    vinv = la.inv(mp.vandermonde(basis.functions, group.unit_nodes))
+    diff = mp.differentiation_matrices(
+            basis.functions, basis.gradients, group.unit_nodes)
     if not isinstance(diff, tuple):
         diff = (diff,)
 
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index ed5fefe53..27ea045ac 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -333,13 +333,11 @@ def make_group_from_vertices(vertices, vertex_indices, order,
 
         # dim, nunit_nodes
         if unit_nodes is None:
-            if dim <= 3:
-                unit_nodes = mp.warp_and_blend_nodes(dim, order)
-            else:
-                unit_nodes = mp.equidistant_nodes(dim, order)
+            shape = mp.Simplex(dim)
+            space = mp.space_for_shape(shape, order)
+            unit_nodes = mp.edge_clustered_nodes_for_space(space, shape)
 
         unit_nodes_01 = 0.5 + 0.5*unit_nodes
-
         nodes = np.einsum(
                 "si,des->dei",
                 unit_nodes_01, spanning_vectors) + el_origins
@@ -352,22 +350,22 @@ def make_group_from_vertices(vertices, vertex_indices, order,
             raise ValueError("invalid number of vertices for tensor-product "
                     "elements, must be power of two")
 
+        shape = mp.Hypercube(dim)
+        space = mp.space_for_shape(shape, order)
+
         if unit_nodes is None:
-            from modepy.quadrature.jacobi_gauss import legendre_gauss_lobatto_nodes
-            unit_nodes = mp.tensor_product_nodes(dim,
-                    legendre_gauss_lobatto_nodes(order))
+            unit_nodes = mp.edge_clustered_nodes_for_space(space, shape)
 
         # shape: (dim, nnodes)
         unit_nodes_01 = 0.5 + 0.5*unit_nodes
         _, nnodes = unit_nodes.shape
 
-        from pytools import generate_nonnegative_integer_tuples_below as gnitb
-        id_tuples = list(gnitb(2, dim))
-        assert len(id_tuples) == nvertices
+        vertex_tuples = mp.node_tuples_for_space(type(space)(dim, 1))
+        assert len(vertex_tuples) == nvertices
 
         vdm = np.empty((nvertices, nvertices))
-        for i, vertex_tuple in enumerate(id_tuples):
-            for j, func_tuple in enumerate(id_tuples):
+        for i, vertex_tuple in enumerate(vertex_tuples):
+            for j, func_tuple in enumerate(vertex_tuples):
                 vertex_ref = np.array(vertex_tuple, dtype=np.float64)
                 vdm[i, j] = np.prod(vertex_ref**func_tuple)
 
@@ -377,7 +375,7 @@ def make_group_from_vertices(vertices, vertex_indices, order,
             coeffs[d] = la.solve(vdm, el_vertices[d].T).T
 
         vdm_nodes = np.zeros((nnodes, nvertices))
-        for j, func_tuple in enumerate(id_tuples):
+        for j, func_tuple in enumerate(vertex_tuples):
             vdm_nodes[:, j] = np.prod(
                     unit_nodes_01 ** np.array(func_tuple).reshape(-1, 1),
                     axis=0)
@@ -767,7 +765,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 
     if dim == 1:
         if mesh_type is not None:
-            raise ValueError("unsupported mesh_type")
+            raise ValueError(f"unsupported mesh type: '{mesh_type}'")
 
         for i in range(shape[0]-1):
             # a--b
@@ -801,7 +799,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
             pass
 
         else:
-            raise ValueError("unsupported mesh_type")
+            raise ValueError(f"unsupported mesh type: '{mesh_type}'")
 
         for i in range(shape[0]-1):
             for j in range(shape[1]-1):
@@ -849,8 +847,8 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
 
                     if is_tp:
                         el_vertices.append(
-                                (a000, a001, a010, a011,
-                                    a100, a101, a110, a111))
+                                (a000, a100, a010, a110,
+                                    a001, a101, a011, a111))
 
                     else:
                         el_vertices.append((a000, a100, a010, a001))
@@ -957,7 +955,8 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
       on [a,b].
     :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
         See :func:`generate_box_mesh`.
-    :param mesh_type: See :func:`generate_box_mesh`.
+    :param group_cls: see :func:`generate_box_mesh`.
+    :param mesh_type: see :func:`generate_box_mesh`.
     """
     if min(n) < 2:
         raise ValueError("need at least two points in each direction")
@@ -975,7 +974,7 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
 
 # {{{ generate_warped_rect_mesh
 
-def generate_warped_rect_mesh(dim, order, n):
+def generate_warped_rect_mesh(dim, order, n, *, group_cls=None):
     """Generate a mesh of a warped line/square/cube. Mainly useful for testing
     functionality with curvilinear meshes.
     """
@@ -983,7 +982,7 @@ def generate_warped_rect_mesh(dim, order, n):
     assert dim in [1, 2, 3]
     mesh = generate_regular_rect_mesh(
             a=(-0.5,)*dim, b=(0.5,)*dim,
-            n=(n,)*dim, order=order)
+            n=(n,)*dim, order=order, group_cls=group_cls)
 
     def m(x):
         result = np.empty_like(x)
diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py
index cf34f407c..d0ac94f7d 100644
--- a/meshmode/mesh/io.py
+++ b/meshmode/mesh/io.py
@@ -191,13 +191,23 @@ def get_mesh(self):
                     continue
 
                 nodes[:, i] = self.points[el_nodes].T
-                vertex_indices[i] = [vertex_gmsh_index_to_mine[v_nr]
-                        for v_nr in el_vertices]
+                vertex_indices[i] = [
+                        vertex_gmsh_index_to_mine[v_nr] for v_nr in el_vertices
+                        ]
 
                 i += 1
 
-            unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(),
-                    dtype=np.float64).T/group_el_type.order)*2 - 1
+            import modepy as mp
+            if isinstance(group_el_type, GmshSimplexElementBase):
+                shape = mp.Simplex(group_el_type.dimensions)
+            elif isinstance(group_el_type, GmshTensorProductElementBase):
+                shape = mp.Hypercube(group_el_type.dimensions)
+            else:
+                raise NotImplementedError(
+                        f"gmsh element type: {type(group_el_type).__name__}")
+
+            space = mp.space_for_shape(shape, group_el_type.order)
+            unit_nodes = mp.equispaced_nodes_for_space(space, shape)
 
             if isinstance(group_el_type, GmshSimplexElementBase):
                 group = SimplexElementGroup(
@@ -213,16 +223,8 @@ def get_mesh(self):
                             np.ones(ngroup_elements, np.bool))
 
             elif isinstance(group_el_type, GmshTensorProductElementBase):
-                gmsh_vertex_tuples = type(group_el_type)(order=1).gmsh_node_tuples()
-                gmsh_vertex_tuples_loc_dict = {
-                        gvt: i
-                        for i, gvt in enumerate(gmsh_vertex_tuples)}
-
-                from pytools import (
-                        generate_nonnegative_integer_tuples_below as gnitb)
-                vertex_shuffle = np.array([
-                    gmsh_vertex_tuples_loc_dict[vt]
-                    for vt in gnitb(2, group_el_type.dimensions)])
+                vertex_shuffle = type(group_el_type)(
+                        order=1).get_lexicographic_gmsh_node_indices()
 
                 group = TensorProductElementGroup(
                     group_el_type.order,
@@ -231,8 +233,8 @@ def get_mesh(self):
                     unit_nodes=unit_nodes
                     )
             else:
-                raise NotImplementedError("gmsh element type: %s"
-                        % type(group_el_type).__name__)
+                # NOTE: already checked above
+                assert False
 
             groups.append(group)
 
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 8644f9326..2054ddd07 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -144,20 +144,20 @@ def Get_size(self):  # noqa: N802
     assert(filecmp.cmp("ref-"+pvtu_filename, pvtu_filename))
 
 
-@pytest.mark.parametrize(("dim", "group_factory"), [
+@pytest.mark.parametrize(("dim", "group_cls"), [
     (1, SimplexElementGroup),
     (2, SimplexElementGroup),
     (3, SimplexElementGroup),
     (2, TensorProductElementGroup),
     (3, TensorProductElementGroup),
     ])
-def test_visualizers(actx_factory, dim, group_factory):
+def test_visualizers(actx_factory, dim, group_cls):
     actx = actx_factory()
 
     nelements = 64
     target_order = 4
 
-    is_simplex = issubclass(group_factory, SimplexElementGroup)
+    is_simplex = issubclass(group_cls, SimplexElementGroup)
     if dim == 1:
         mesh = mgen.make_curve_mesh(
                 mgen.NArmedStarfish(5, 0.25),
@@ -169,7 +169,7 @@ def test_visualizers(actx_factory, dim, group_factory):
         else:
             mesh = mgen.generate_regular_rect_mesh(
                     a=(0,)*dim, b=(1,)*dim, n=(5,)*dim,
-                    group_factory=group_factory,
+                    group_cls=group_cls,
                     order=target_order)
     elif dim == 3:
         if is_simplex:
@@ -177,18 +177,18 @@ def test_visualizers(actx_factory, dim, group_factory):
         else:
             mesh = mgen.generate_regular_rect_mesh(
                     a=(0,)*dim, b=(1,)*dim, n=(5,)*dim,
-                    group_factory=group_factory,
+                    group_cls=group_cls,
                     order=target_order)
     else:
         raise ValueError("unknown dimensionality")
 
     if is_simplex:
-        discr_group_factory = InterpolatoryQuadratureSimplexGroupFactory
+        group_factory = InterpolatoryQuadratureSimplexGroupFactory
     else:
-        discr_group_factory = LegendreGaussLobattoTensorProductGroupFactory
+        group_factory = LegendreGaussLobattoTensorProductGroupFactory
 
     from meshmode.discretization import Discretization
-    discr = Discretization(actx, mesh, discr_group_factory(target_order))
+    discr = Discretization(actx, mesh, group_factory(target_order))
 
     nodes = thaw(actx, discr.nodes())
     f = actx.np.sqrt(sum(nodes**2)) + 1j*nodes[0]
@@ -285,14 +285,16 @@ def test_boundary_tags():
     ])
 @pytest.mark.parametrize("group_cls", [
     SimplexElementGroup,
-
-    # FIXME: Not implemented: TPE.face_vertex_indices
-    # TensorProductElementGroup
+    TensorProductElementGroup
     ])
 def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False):
+    if group_cls is TensorProductElementGroup and mesh_type is not None:
+        pytest.skip("mesh type not supported on tensor product elements")
+
     from meshmode.mesh.generation import generate_regular_rect_mesh
     from meshmode.mesh import is_boundary_tag_empty
     from meshmode.mesh import check_bc_coverage
+
     if dim == 1:
         a = (0,)
         b = (1,)
@@ -317,7 +319,7 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False):
                                       group_cls=group_cls,
                                       mesh_type=mesh_type)
 
-    if visualize:
+    if visualize and dim == 2:
         from meshmode.mesh.visualization import draw_2d_mesh
         draw_2d_mesh(mesh, draw_element_numbers=False, draw_vertex_numbers=False)
         import matplotlib.pyplot as plt
@@ -326,8 +328,12 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False):
     # correct answer
     if dim == 1:
         num_on_bdy = 1
+    elif group_cls is TensorProductElementGroup:
+        num_on_bdy = dim * (nelem-1)**(dim-1)
+    elif group_cls is SimplexElementGroup:
+        num_on_bdy = dim * (dim-1) * (nelem-1)**(dim-1)
     else:
-        num_on_bdy = dim*(dim-1)*(nelem-1)**(dim-1)
+        assert False
 
     assert not is_boundary_tag_empty(mesh, "btag_test_1")
     assert not is_boundary_tag_empty(mesh, "btag_test_2")
@@ -358,7 +364,6 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False):
         raise ValueError("%i marked on custom boundary 2, should be %i" %
                          (num_marked_bdy_2, num_on_bdy))
 
-
 # }}}
 
 
@@ -368,7 +373,11 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False):
     InterpolatoryQuadratureSimplexGroupFactory,
     PolynomialWarpAndBlendGroupFactory,
     partial(PolynomialRecursiveNodesGroupFactory, family="lgl"),
-    #partial(PolynomialRecursiveNodesGroupFactory, family="gc"),
+
+    # Redundant, no information gain.
+    # partial(PolynomialRecursiveNodesGroupFactory, family="gc"),
+
+    LegendreGaussLobattoTensorProductGroupFactory,
     ])
 @pytest.mark.parametrize("boundary_tag", [
     BTAG_ALL,
@@ -377,14 +386,27 @@ def test_box_boundary_tags(dim, nelem, mesh_type, group_cls, visualize=False):
     ])
 @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
     ("blob", 2, ["8e-2", "6e-2", "4e-2"]),
+
+    # If "warp" works, "rect" likely does, too.
+    # ("rect", 3, [10, 20, 30]),
+
     ("warp", 2, [10, 20, 30]),
     ("warp", 3, [10, 20, 30]),
     ])
-@pytest.mark.parametrize("per_face_groups", [False, True])
+@pytest.mark.parametrize("per_face_groups", [True, False])
 def test_boundary_interpolation(actx_factory, group_factory, boundary_tag,
         mesh_name, dim, mesh_pars, per_face_groups):
+    if (group_factory is LegendreGaussLobattoTensorProductGroupFactory
+            and mesh_name == "blob"):
+        pytest.skip("tensor products not implemented on blobs")
+
     actx = actx_factory()
 
+    if group_factory is LegendreGaussLobattoTensorProductGroupFactory:
+        group_cls = TensorProductElementGroup
+    else:
+        group_cls = SimplexElementGroup
+
     from meshmode.discretization import Discretization
     from meshmode.discretization.connection import (
             make_face_restriction, check_connection)
@@ -420,7 +442,15 @@ def f(x):
                     force_ambient_dim=2)
         elif mesh_name == "warp":
             from meshmode.mesh.generation import generate_warped_rect_mesh
-            mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par)
+            mesh = generate_warped_rect_mesh(dim, order=order, n=mesh_par,
+                    group_cls=group_cls)
+
+            h = 1/mesh_par
+
+        elif mesh_name == "rect":
+            from meshmode.mesh.generation import generate_regular_rect_mesh
+            mesh = generate_regular_rect_mesh(a=(0,)*dim, b=(1,)*dim,
+                    order=order, n=(mesh_par,)*dim, group_cls=group_cls)
 
             h = 1/mesh_par
         else:
@@ -462,23 +492,37 @@ def f(x):
     print(eoc_rec)
     assert (
             eoc_rec.order_estimate() >= order-order_slack
-            or eoc_rec.max_error() < 3e-14)
+            or eoc_rec.max_error() < 1.5e-13)
 
 # }}}
 
 
 # {{{ boundary-to-all-faces connecttion
 
+@pytest.mark.parametrize("group_factory", [
+    InterpolatoryQuadratureSimplexGroupFactory,
+    PolynomialWarpAndBlendGroupFactory,
+    LegendreGaussLobattoTensorProductGroupFactory,
+    ])
 @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
     ("blob", 2, [1e-1, 8e-2, 5e-2]),
     ("warp", 2, [10, 20, 30]),
     ("warp", 3, [10, 20, 30]),
     ])
 @pytest.mark.parametrize("per_face_groups", [False, True])
-def test_all_faces_interpolation(actx_factory, mesh_name, dim, mesh_pars,
-        per_face_groups):
+def test_all_faces_interpolation(actx_factory, group_factory,
+        mesh_name, dim, mesh_pars, per_face_groups):
+    if (group_factory is LegendreGaussLobattoTensorProductGroupFactory
+            and mesh_name == "blob"):
+        pytest.skip("tensor products not implemented on blobs")
+
     actx = actx_factory()
 
+    if group_factory is LegendreGaussLobattoTensorProductGroupFactory:
+        group_cls = TensorProductElementGroup
+    else:
+        group_cls = SimplexElementGroup
+
     from meshmode.discretization import Discretization
     from meshmode.discretization.connection import (
             make_face_restriction, make_face_to_all_faces_embedding,
@@ -512,7 +556,8 @@ def f(x):
             print("END GEN")
         elif mesh_name == "warp":
             from meshmode.mesh.generation import generate_warped_rect_mesh
-            mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par)
+            mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par,
+                    group_cls=group_cls)
 
             h = 1/mesh_par
         else:
@@ -520,13 +565,12 @@ def f(x):
 
         # }}}
 
-        vol_discr = Discretization(actx, mesh,
-                PolynomialWarpAndBlendGroupFactory(order))
+        vol_discr = Discretization(actx, mesh, group_factory(order))
         print("h=%s -> %d elements" % (
                 h, sum(mgrp.nelements for mgrp in mesh.groups)))
 
         all_face_bdry_connection = make_face_restriction(
-                actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order),
+                actx, vol_discr, group_factory(order),
                 FACE_RESTR_ALL, per_face_groups=per_face_groups)
         all_face_bdry_discr = all_face_bdry_connection.to_discr
 
@@ -551,7 +595,7 @@ def f(x):
                 FACE_RESTR_INTERIOR,
                 ]:
             bdry_connection = make_face_restriction(
-                    actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order),
+                    actx, vol_discr, group_factory(order),
                     boundary_tag, per_face_groups=per_face_groups)
             bdry_discr = bdry_connection.to_discr
 
@@ -580,19 +624,29 @@ def f(x):
 
 @pytest.mark.parametrize("group_factory", [
     InterpolatoryQuadratureSimplexGroupFactory,
-    PolynomialWarpAndBlendGroupFactory
+    PolynomialWarpAndBlendGroupFactory,
+    LegendreGaussLobattoTensorProductGroupFactory,
     ])
 @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [
     ("segment", 1, [8, 16, 32]),
     ("blob", 2, [1e-1, 8e-2, 5e-2]),
     ("warp", 2, [3, 5, 7]),
-    ("warp", 3, [3, 5]),
+    ("warp", 3, [5, 7])
     ])
 def test_opposite_face_interpolation(actx_factory, group_factory,
         mesh_name, dim, mesh_pars):
+    if (group_factory is LegendreGaussLobattoTensorProductGroupFactory
+            and mesh_name in ["segment", "blob"]):
+        pytest.skip("tensor products not implemented on blobs")
+
     logging.basicConfig(level=logging.INFO)
     actx = actx_factory()
 
+    if group_factory is LegendreGaussLobattoTensorProductGroupFactory:
+        group_cls = TensorProductElementGroup
+    else:
+        group_cls = SimplexElementGroup
+
     from meshmode.discretization import Discretization
     from meshmode.discretization.connection import (
             make_face_restriction, make_opposite_face_connection,
@@ -615,7 +669,8 @@ def f(x):
             from meshmode.mesh.generation import generate_box_mesh
             mesh = generate_box_mesh(
                     [np.linspace(-0.5, 0.5, mesh_par)],
-                    order=order)
+                    order=order,
+                    group_cls=group_cls)
             h = 1.0 / mesh_par
         elif mesh_name == "blob":
             assert dim == 2
@@ -634,7 +689,8 @@ def f(x):
             print("END GEN")
         elif mesh_name == "warp":
             from meshmode.mesh.generation import generate_warped_rect_mesh
-            mesh = generate_warped_rect_mesh(dim, order=4, n=mesh_par)
+            mesh = generate_warped_rect_mesh(dim, order=order, n=mesh_par,
+                    group_cls=group_cls)
 
             h = 1/mesh_par
         else:
@@ -642,8 +698,7 @@ def f(x):
 
         # }}}
 
-        vol_discr = Discretization(actx, mesh,
-                group_factory(order))
+        vol_discr = Discretization(actx, mesh, group_factory(order))
         print("h=%s -> %d elements" % (
                 h, sum(mgrp.nelements for mgrp in mesh.groups)))
 
@@ -665,7 +720,7 @@ def f(x):
     print(eoc_rec)
     assert (
             eoc_rec.order_estimate() >= order-0.5
-            or eoc_rec.max_error() < 1e-13)
+            or eoc_rec.max_error() < 1.7e-13)
 
 # }}}
 
@@ -767,9 +822,7 @@ def test_orientation_3d(actx_factory, what, mesh_gen_func, visualize=False):
 
 @pytest.mark.parametrize("group_cls", [
     SimplexElementGroup,
-
-    # NOTE: needs TensorProductElementGroup.face_vertex_indices
-    # TensorProductElementGroup
+    TensorProductElementGroup
     ])
 def test_merge_and_map(actx_factory, group_cls, visualize=False):
     from meshmode.mesh.io import generate_gmsh, FileSource
@@ -822,43 +875,57 @@ def test_merge_and_map(actx_factory, group_cls, visualize=False):
 # {{{ sanity checks: single element
 
 @pytest.mark.parametrize("dim", [2, 3])
-@pytest.mark.parametrize("order", [1, 3])
-def test_sanity_single_element(actx_factory, dim, order, visualize=False):
+@pytest.mark.parametrize("mesh_order", [1, 3])
+@pytest.mark.parametrize("group_cls", [
+    SimplexElementGroup,
+    TensorProductElementGroup,
+    ])
+def test_sanity_single_element(actx_factory, dim, mesh_order, group_cls,
+        visualize=False):
     pytest.importorskip("pytential")
     actx = actx_factory()
 
-    from modepy.tools import unit_vertices
-    vertices = unit_vertices(dim).T.copy()
+    if group_cls is SimplexElementGroup:
+        from meshmode.discretization.poly_element import \
+                PolynomialWarpAndBlendGroupFactory as GroupFactory
+    elif group_cls is TensorProductElementGroup:
+        from meshmode.discretization.poly_element import \
+                LegendreGaussLobattoTensorProductGroupFactory as GroupFactory
+    else:
+        raise TypeError
+
+    import modepy as mp
+    shape = group_cls._modepy_shape_cls(dim)
+    space = mp.space_for_shape(shape, mesh_order)
+
+    vertices = mp.unit_vertices_for_shape(shape)
+    nodes = mp.edge_clustered_nodes_for_space(space, shape).reshape(dim, 1, -1)
+    vertex_indices = np.arange(shape.nvertices, dtype=np.int32).reshape(1, -1)
 
     center = np.empty(dim, np.float64)
     center.fill(-0.5)
 
-    import modepy as mp
-    from meshmode.mesh import SimplexElementGroup, Mesh, BTAG_ALL
-    mg = SimplexElementGroup(
-            order=order,
-            vertex_indices=np.arange(dim+1, dtype=np.int32).reshape(1, -1),
-            nodes=mp.warp_and_blend_nodes(dim, order).reshape(dim, 1, -1),
-            dim=dim)
-
+    from meshmode.mesh import Mesh, BTAG_ALL
+    mg = group_cls(mesh_order, vertex_indices, nodes, dim=dim)
     mesh = Mesh(vertices, [mg], is_conforming=True)
 
     from meshmode.discretization import Discretization
-    from meshmode.discretization.poly_element import \
-            PolynomialWarpAndBlendGroupFactory
-    vol_discr = Discretization(actx, mesh,
-            PolynomialWarpAndBlendGroupFactory(order+3))
+    vol_discr = Discretization(actx, mesh, GroupFactory(mesh_order + 3))
 
     # {{{ volume calculation check
 
-    vol_x = thaw(actx, vol_discr.nodes())
-
-    vol_one = vol_x[0] * 0 + 1
-    from pytential import norm, integral  # noqa
+    if isinstance(mg, SimplexElementGroup):
+        from pytools import factorial
+        true_vol = 1/factorial(dim) * 2**dim
+    elif isinstance(mg, TensorProductElementGroup):
+        true_vol = 2**dim
+    else:
+        raise TypeError
 
-    from pytools import factorial
-    true_vol = 1/factorial(dim) * 2**dim
+    nodes = thaw(actx, vol_discr.nodes())
+    vol_one = 1 + 0 * nodes[0]
 
+    from pytential import norm, integral  # noqa
     comp_vol = integral(vol_discr, vol_one)
     rel_vol_err = abs(true_vol - comp_vol) / true_vol
 
@@ -870,14 +937,14 @@ def test_sanity_single_element(actx_factory, dim, order, visualize=False):
 
     from meshmode.discretization.connection import make_face_restriction
     bdry_connection = make_face_restriction(
-            actx, vol_discr, PolynomialWarpAndBlendGroupFactory(order + 3),
+            actx, vol_discr, GroupFactory(mesh_order + 3),
             BTAG_ALL)
     bdry_discr = bdry_connection.to_discr
 
     # }}}
 
     from pytential import bind, sym
-    bdry_normals = bind(bdry_discr, sym.normal(dim))(actx).as_vector(dtype=object)
+    bdry_normals = bind(bdry_discr, sym.normal(dim).as_vector())(actx)
 
     if visualize:
         from meshmode.discretization.visualization import make_visualizer
@@ -1210,39 +1277,70 @@ def test_quad_mesh_2d(ambient_dim, filename, visualize=False):
 
 # {{{ test_quad_mesh_3d
 
-# This currently (gmsh 2.13.2) crashes gmsh. A massaged version of this using
-# 'cube.step' succeeded in generating 'hybrid-cube.msh' and 'cubed-cube.msh'.
-def no_test_quad_mesh_3d():
-    from meshmode.mesh.io import generate_gmsh, ScriptWithFilesSource
-    print("BEGIN GEN")
-    mesh = generate_gmsh(
-            ScriptWithFilesSource(
-                """
-                Merge "ball-radius-1.step";
-                // Mesh.CharacteristicLengthMax = 0.1;
+@pytest.mark.parametrize("mesh_name", [
+    # this currently (2020-11-05 with gmsh 4.6.0) does not recombine anything
+    # or flat out crashes gmsh
+    # "ball",
 
-                Mesh.RecombineAll=1;
-                Mesh.Recombine3DAll=1;
-                Mesh.Algorithm = 8;
-                Mesh.Algorithm3D = 9;
-                // Mesh.Smoothing = 0;
+    "cube",
+    ])
+def test_quad_mesh_3d(mesh_name, order=3, visualize=False):
+    if mesh_name == "ball":
+        from meshmode.mesh.io import ScriptWithFilesSource
+        script = ScriptWithFilesSource(
+            """
+            Merge "ball-radius-1.step";
+            // Mesh.CharacteristicLengthMax = 0.1;
+
+            Mesh.RecombineAll = 1;
+            Mesh.Recombine3DAll = 1;
+            Mesh.Recombine3DLevel = 2;
+
+            Mesh.Algorithm = 8;
+            Mesh.Algorithm3D = 8;
+
+            Mesh 3;
+            Save "output.msh";
+            """,
+            ["ball-radius-1.step"])
+
+        with open("ball-quad.geo", "w") as f:
+            f.write(script.source)
+
+    elif mesh_name == "cube":
+        from meshmode.mesh.io import ScriptSource
+        script = ScriptSource(
+            """
+            SetFactory("OpenCASCADE");
+            Box(1) = {0, 0, 0, 1, 1, 1};
+
+            Transfinite Line "*" = 8;
+            Transfinite Surface "*";
+            Transfinite Volume "*";
+
+            Mesh.RecombineAll = 1;
+            Mesh.Recombine3DAll = 1;
+            Mesh.Recombine3DLevel = 2;
+            """, "geo")
+    else:
+        raise ValueError(f"unknown mesh name: '{mesh_name}'")
 
-                // Mesh.ElementOrder = 3;
+    np.set_printoptions(linewidth=200)
+    from meshmode.mesh.io import generate_gmsh
+    logger.info("BEGIN GEN")
+    mesh = generate_gmsh(script, 3, order=order, target_unit="MM")
+    logger.info("END GEN")
 
-                Mesh 3;
-                Save "output.msh";
-                """,
-                ["ball-radius-1.step"]),
-            )
-    print("END GEN")
-    print(mesh.nelements)
+    if visualize:
+        from meshmode.mesh.visualization import write_vertex_vtk_file
+        write_vertex_vtk_file(mesh, f"quad_mesh_3d_{mesh_name}.vtu", overwrite=True)
 
 # }}}
 
 
 # {{{ test_quad_single_element
 
-def test_quad_single_element():
+def test_quad_single_element(visualize=False):
     from meshmode.mesh.generation import make_group_from_vertices
 
     vertices = np.array([
@@ -1257,7 +1355,7 @@ def test_quad_single_element():
             30, group_cls=TensorProductElementGroup)
 
     Mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None)
-    if 0:
+    if visualize:
         import matplotlib.pyplot as plt
         plt.plot(
                 mg.nodes[0].reshape(-1),