diff --git a/FIAT/__init__.py b/FIAT/__init__.py
index 3c99570e0..26c467a99 100644
--- a/FIAT/__init__.py
+++ b/FIAT/__init__.py
@@ -48,9 +48,12 @@
 from FIAT.raviart_thomas import RaviartThomas
 from FIAT.crouzeix_raviart import CrouzeixRaviart
 from FIAT.regge import Regge
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson
 from FIAT.arnold_winther import ArnoldWinther
 from FIAT.arnold_winther import ArnoldWintherNC
+from FIAT.hu_zhang import HuZhang
 from FIAT.mardal_tai_winther import MardalTaiWinther
 from FIAT.bubble import Bubble, FacetBubble
 from FIAT.tensor_product import TensorProductElement
@@ -107,8 +110,11 @@
                       "BrokenElement": DiscontinuousElement,
                       "HDiv Trace": HDivTrace,
                       "Hellan-Herrmann-Johnson": HellanHerrmannJohnson,
+                      "Gopalakrishnan-Lederer-Schoberl 1st kind": GopalakrishnanLedererSchoberlFirstKind,
+                      "Gopalakrishnan-Lederer-Schoberl 2nd kind": GopalakrishnanLedererSchoberlSecondKind,
                       "Conforming Arnold-Winther": ArnoldWinther,
                       "Nonconforming Arnold-Winther": ArnoldWintherNC,
+                      "Hu-Zhang": HuZhang,
                       "Mardal-Tai-Winther": MardalTaiWinther}
 
 # List of extra elements
diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py
index 5c12fcc76..e7af81783 100644
--- a/FIAT/arnold_winther.py
+++ b/FIAT/arnold_winther.py
@@ -9,157 +9,125 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 
-from FIAT.finite_element import CiarletElement
-from FIAT.dual_set import DualSet
-from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet
-from FIAT.functional import (
-    PointwiseInnerProductEvaluation as InnerProduct,
-    FrobeniusIntegralMoment as FIM,
-    IntegralMomentOfTensorDivergence,
-    IntegralLegendreNormalNormalMoment,
-    IntegralLegendreNormalTangentialMoment)
-
-from FIAT.quadrature import make_quadrature
+from FIAT import finite_element, dual_set, polynomial_set
+from FIAT.reference_element import TRIANGLE
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.functional import (ComponentPointEvaluation,
+                             TensorBidirectionalIntegralMoment,
+                             IntegralMomentOfTensorDivergence,
+                             IntegralLegendreNormalNormalMoment,
+                             IntegralLegendreNormalTangentialMoment)
 
 import numpy
 
 
-class ArnoldWintherNCDual(DualSet):
-    def __init__(self, cell, degree=2):
-        if not degree == 2:
-            raise ValueError("Nonconforming Arnold-Winther elements are"
-                             "only defined for degree 2.")
-        dofs = []
-        dof_ids = {}
-        dof_ids[0] = {0: [], 1: [], 2: []}
-        dof_ids[1] = {0: [], 1: [], 2: []}
-        dof_ids[2] = {0: []}
+class ArnoldWintherNCDual(dual_set.DualSet):
+    def __init__(self, ref_el, degree=2):
+        if degree != 2:
+            raise ValueError("Nonconforming Arnold-Winther elements are only defined for degree 2.")
+        top = ref_el.get_topology()
+        sd = ref_el.get_spatial_dimension()
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+        nodes = []
 
-        dof_cur = 0
         # no vertex dofs
         # proper edge dofs now (not the contraints)
-        # moments of normal . sigma against constants and linears.
-        for entity_id in range(3):                  # a triangle has 3 edges
-            for order in (0, 1):
-                dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6),
-                         IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)]
-            dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4))
-            dof_cur += 4
+        # edge dofs: bidirectional nn and nt moments against P1.
+        qdegree = degree + 2
+        for entity in sorted(top[1]):
+            cur = len(nodes)
+            for order in range(2):
+                nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree))
+                nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree))
+            entity_ids[1][entity].extend(range(cur, len(nodes)))
 
         # internal dofs: constant moments of three unique components
-        Q = make_quadrature(cell, 2)
-
-        e1 = numpy.array([1.0, 0.0])              # euclidean basis 1
-        e2 = numpy.array([0.0, 1.0])              # euclidean basis 2
-        basis = [(e1, e1), (e1, e2), (e2, e2)]    # basis for symmetric matrices
-        for (v1, v2) in basis:
-            v1v2t = numpy.outer(v1, v2)
-            fatqp = numpy.zeros((2, 2, len(Q.pts)))
-            for i, y in enumerate(v1v2t):
-                for j, x in enumerate(y):
-                    for k in range(len(Q.pts)):
-                        fatqp[i, j, k] = x
-            dofs.append(FIM(cell, Q, fatqp))
-        dof_ids[2][0] = list(range(dof_cur, dof_cur + 3))
-        dof_cur += 3
+        cur = len(nodes)
+        n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
+        Q = create_quadrature(ref_el, degree)
+        phi = numpy.full(Q.get_weights().shape, 1/ref_el.volume())
+        nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
+                     for i in range(sd) for j in range(i, sd))
+        entity_ids[2][0].extend(range(cur, len(nodes)))
 
         # put the constraint dofs last.
-        for entity_id in range(3):
-            dof = IntegralLegendreNormalNormalMoment(cell, entity_id, 2, 6)
-            dofs.append(dof)
-            dof_ids[1][entity_id].append(dof_cur)
-            dof_cur += 1
+        for entity in sorted(top[1]):
+            cur = len(nodes)
+            nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, 2, qdegree))
+            entity_ids[1][entity].append(cur)
 
-        super().__init__(dofs, cell, dof_ids)
+        super().__init__(nodes, ref_el, entity_ids)
 
 
-class ArnoldWintherNC(CiarletElement):
+class ArnoldWintherNC(finite_element.CiarletElement):
     """The definition of the nonconforming Arnold-Winther element.
     """
-    def __init__(self, cell, degree=2):
-        assert degree == 2, "Only defined for degree 2"
-        Ps = ONSymTensorPolynomialSet(cell, degree)
-        Ls = ArnoldWintherNCDual(cell, degree)
+    def __init__(self, ref_el, degree=2):
+        if ref_el.shape != TRIANGLE:
+            raise ValueError(f"{type(self).__name__} only defined on triangles")
+        Ps = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree)
+        Ls = ArnoldWintherNCDual(ref_el, degree)
+        formdegree = ref_el.get_spatial_dimension() - 1
         mapping = "double contravariant piola"
+        super().__init__(Ps, Ls, degree, formdegree, mapping=mapping)
 
-        super().__init__(Ps, Ls, degree, mapping=mapping)
 
-
-class ArnoldWintherDual(DualSet):
-    def __init__(self, cell, degree=3):
-        if not degree == 3:
-            raise ValueError("Arnold-Winther elements are"
-                             "only defined for degree 3.")
-        dofs = []
-        dof_ids = {}
-        dof_ids[0] = {0: [], 1: [], 2: []}
-        dof_ids[1] = {0: [], 1: [], 2: []}
-        dof_ids[2] = {0: []}
-
-        dof_cur = 0
+class ArnoldWintherDual(dual_set.DualSet):
+    def __init__(self, ref_el, degree=3):
+        if degree != 3:
+            raise ValueError("Arnold-Winther elements are only defined for degree 3.")
+        top = ref_el.get_topology()
+        sd = ref_el.get_spatial_dimension()
+        shp = (sd, sd)
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+        nodes = []
 
         # vertex dofs
-        vs = cell.get_vertices()
-        e1 = numpy.array([1.0, 0.0])
-        e2 = numpy.array([0.0, 1.0])
-        basis = [(e1, e1), (e1, e2), (e2, e2)]
-
-        dof_cur = 0
-
-        for entity_id in range(3):
-            node = tuple(vs[entity_id])
-            for (v1, v2) in basis:
-                dofs.append(InnerProduct(cell, v1, v2, node))
-            dof_ids[0][entity_id] = list(range(dof_cur, dof_cur + 3))
-            dof_cur += 3
-
-        # edge dofs now
-        # moments of normal . sigma against constants and linears.
-        for entity_id in range(3):
-            for order in (0, 1):
-                dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6),
-                         IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)]
-            dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4))
-            dof_cur += 4
-
-        # internal dofs: constant moments of three unique components
-        Q = make_quadrature(cell, 3)
-
-        e1 = numpy.array([1.0, 0.0])              # euclidean basis 1
-        e2 = numpy.array([0.0, 1.0])              # euclidean basis 2
-        basis = [(e1, e1), (e1, e2), (e2, e2)]    # basis for symmetric matrices
-        for (v1, v2) in basis:
-            v1v2t = numpy.outer(v1, v2)
-            fatqp = numpy.zeros((2, 2, len(Q.pts)))
-            for k in range(len(Q.pts)):
-                fatqp[:, :, k] = v1v2t
-            dofs.append(FIM(cell, Q, fatqp))
-        dof_ids[2][0] = list(range(dof_cur, dof_cur + 3))
-        dof_cur += 3
-
-        # Constraint dofs
-
-        Q = make_quadrature(cell, 5)
-
-        onp = ONPolynomialSet(cell, 2, (2,))
-        pts = Q.get_points()
-        onpvals = onp.tabulate(pts)[0, 0]
-
-        for i in list(range(3, 6)) + list(range(9, 12)):
-            dofs.append(IntegralMomentOfTensorDivergence(cell, Q,
-                                                         onpvals[i, :, :]))
-
-        dof_ids[2][0] += list(range(dof_cur, dof_cur+6))
-
-        super().__init__(dofs, cell, dof_ids)
-
-
-class ArnoldWinther(CiarletElement):
+        for v in sorted(top[0]):
+            cur = len(nodes)
+            pt, = ref_el.make_points(0, v, degree)
+            nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt)
+                         for i in range(sd) for j in range(i, sd))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        # edge dofs: bidirectional nn and nt moments against P_{k-2}
+        max_order = degree - 2
+        qdegree = degree + max_order
+        for entity in sorted(top[1]):
+            cur = len(nodes)
+            for order in range(max_order+1):
+                nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree))
+                nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree))
+            entity_ids[1][entity].extend(range(cur, len(nodes)))
+
+        # internal dofs: moments of unique components against P_{k-3}
+        n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
+        Q = create_quadrature(ref_el, 2*(degree-1))
+        P = polynomial_set.ONPolynomialSet(ref_el, degree-3, scale="L2 piola")
+        phis = P.tabulate(Q.get_points())[(0,)*sd]
+        nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
+                     for phi in phis for i in range(sd) for j in range(i, sd))
+
+        # constraint dofs: moments of divergence against P_{k-1} \ P_{k-2}
+        P = polynomial_set.ONPolynomialSet(ref_el, degree-1, shape=(sd,))
+        dimPkm1 = P.expansion_set.get_num_members(degree-1)
+        dimPkm2 = P.expansion_set.get_num_members(degree-2)
+        PH = P.take([i + j * dimPkm1 for j in range(sd) for i in range(dimPkm2, dimPkm1)])
+        phis = PH.tabulate(Q.get_points())[(0,)*sd]
+        nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis)
+
+        entity_ids[2][0].extend(range(cur, len(nodes)))
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class ArnoldWinther(finite_element.CiarletElement):
     """The definition of the conforming Arnold-Winther element.
     """
-    def __init__(self, cell, degree=3):
-        assert degree == 3, "Only defined for degree 3"
-        Ps = ONSymTensorPolynomialSet(cell, degree)
-        Ls = ArnoldWintherDual(cell, degree)
+    def __init__(self, ref_el, degree=3):
+        if ref_el.shape != TRIANGLE:
+            raise ValueError(f"{type(self).__name__} only defined on triangles")
+        Ps = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree)
+        Ls = ArnoldWintherDual(ref_el, degree)
+        formdegree = ref_el.get_spatial_dimension() - 1
         mapping = "double contravariant piola"
-        super().__init__(Ps, Ls, degree, mapping=mapping)
+        super().__init__(Ps, Ls, degree, formdegree, mapping=mapping)
diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index 7fa0ad887..28eb561d4 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -279,16 +279,19 @@ def __init__(self, ref_el, scale=None, variant=None):
         self._dmats_cache = {}
         self._cell_node_map_cache = {}
 
-    def get_scale(self, cell=0):
+    def get_scale(self, n, cell=0):
         scale = self.scale
+        sd = self.ref_el.get_spatial_dimension()
         if isinstance(scale, str):
-            sd = self.ref_el.get_spatial_dimension()
             vol = self.ref_el.volume_of_subcomplex(sd, cell)
             scale = scale.lower()
             if scale == "orthonormal":
                 scale = math.sqrt(1.0 / vol)
             elif scale == "l2 piola":
                 scale = 1.0 / vol
+        elif n == 0 and sd > 1 and len(self.affine_mappings) == 1:
+            # return 1 for n=0 to make regression tests pass
+            scale = 1
         return scale
 
     def get_num_members(self, n):
@@ -310,9 +313,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         ref_pts = numpy.add(numpy.dot(pts, A.T), b).T
         Jinv = A if direction is None else numpy.dot(A, direction)[:, None]
         sd = self.ref_el.get_spatial_dimension()
-
-        # Always return 1 for n=0 to make regression tests pass
-        scale = 1.0 if n == 0 and len(self.affine_mappings) == 1 else self.get_scale(cell=cell)
+        scale = self.get_scale(n, cell=cell)
         phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
                                  scale, variant=self.variant)
         if self.continuity == "C0":
@@ -549,7 +550,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         Jinv = A[0, 0] if direction is None else numpy.dot(A, direction)
         xs = numpy.add(numpy.dot(pts, A.T), b)
         results = {}
-        scale = self.get_scale(cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
+        scale = self.get_scale(n, cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1)
         for k in range(order+1):
             v = numpy.zeros((n + 1, len(xs)), xs.dtype)
             if n >= k:
diff --git a/FIAT/functional.py b/FIAT/functional.py
index 720412c5d..6a831a87c 100644
--- a/FIAT/functional.py
+++ b/FIAT/functional.py
@@ -631,13 +631,13 @@ class TensorBidirectionalIntegralMoment(FrobeniusIntegralMoment):
     field, f a function tabulated at points, and v,w be vectors. This implements the evaluation
     \int v^T u(x) w f(x).
     Clearly v^iu_{ij}w^j = u_{ij}v^iw^j. Thus the value can be computed
-    from the Frobenius inner product of u with wv^T. This gives the
+    from the Frobenius inner product of u with vw^T. This gives the
     correct weights.
     """
 
     def __init__(self, ref_el, v, w, Q, f_at_qpts):
-        wvT = numpy.outer(w, v)
-        F_at_qpts = numpy.multiply(wvT[..., None], f_at_qpts)
+        vwT = numpy.outer(v, w)
+        F_at_qpts = numpy.multiply(vwT[..., None], f_at_qpts)
         super().__init__(ref_el, Q, F_at_qpts, "TensorBidirectionalMomentInnerProductEvaluation")
 
 
diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py
new file mode 100644
index 000000000..656cee6bb
--- /dev/null
+++ b/FIAT/gopalakrishnan_lederer_schoberl.py
@@ -0,0 +1,88 @@
+from FIAT import finite_element, dual_set, polynomial_set, expansions
+from FIAT.functional import TensorBidirectionalIntegralMoment as BidirectionalMoment
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.restricted import RestrictedElement
+
+
+class GLSDual(dual_set.DualSet):
+    def __init__(self, ref_el, degree):
+        sd = ref_el.get_spatial_dimension()
+        top = ref_el.get_topology()
+        nodes = []
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+
+        # Face dofs: moments of normal-tangential components against a basis for Pk
+        ref_facet = ref_el.construct_subelement(sd-1)
+        Qref = create_quadrature(ref_facet, 2*degree)
+        P = polynomial_set.ONPolynomialSet(ref_facet, degree)
+        phis = P.tabulate(Qref.get_points())[(0,) * (sd-1)]
+        for f in sorted(top[sd-1]):
+            cur = len(nodes)
+            Q = FacetQuadratureRule(ref_el, sd-1, f, Qref)
+            Jdet = Q.jacobian_determinant()
+            normal = ref_el.compute_scaled_normal(f)
+            tangents = ref_el.compute_tangents(sd-1, f)
+            n = normal / Jdet
+            nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
+                         for phi in phis for t in tangents)
+            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+
+        # Interior dofs: moments of normal-tangential components against a basis for P_{k-1}
+        if degree > 0:
+            cur = len(nodes)
+            Q = create_quadrature(ref_el, 2*degree-1)
+            P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola")
+            phis = P.tabulate(Q.get_points())[(0,) * sd]
+            for f in sorted(top[sd-1]):
+                n = ref_el.compute_scaled_normal(f)
+                tangents = ref_el.compute_tangents(sd-1, f)
+                nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
+                             for phi in phis for t in tangents)
+            entity_ids[sd][0].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement):
+    """The GLS element used for the Mass-Conserving mixed Stress (MCS)
+    formulation for Stokes flow with weakly imposed stress symmetry.
+
+    GLS^2(k) is the space of trace-free polynomials of degree k with
+    continuous normal-tangential components.
+
+    Reference: https://doi.org/10.1137/19M1248960
+
+    Notes
+    -----
+    This element does not include the bubbles required for inf-sup stability of
+    the weak symmetry constraint.
+
+    """
+    def __init__(self, ref_el, degree):
+        poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree)
+        dual = GLSDual(ref_el, degree)
+        sd = ref_el.get_spatial_dimension()
+        formdegree = (1, sd-1)
+        mapping = "covariant contravariant piola"
+        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
+
+
+def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree):
+    """The GLS element used for the Mass-Conserving mixed Stress (MCS)
+    formulation for Stokes flow.
+
+    GLS^1(k) is the space of trace-free polynomials of degree k with
+    continuous normal-tangential components of degree k-1.
+
+    Reference: https://doi.org/10.1093/imanum/drz022
+    """
+    fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree)
+    entity_dofs = fe.entity_dofs()
+    sd = ref_el.get_spatial_dimension()
+    dimPkm1 = (sd-1)*expansions.polynomial_dimension(ref_el.construct_subelement(sd-1), degree-1)
+    indices = []
+    for f in entity_dofs[sd-1]:
+        indices.extend(entity_dofs[sd-1][f][:dimPkm1])
+    indices.extend(entity_dofs[sd][0])
+    return RestrictedElement(fe, indices=indices)
diff --git a/FIAT/guzman_neilan.py b/FIAT/guzman_neilan.py
index ff853e583..970dbdf89 100644
--- a/FIAT/guzman_neilan.py
+++ b/FIAT/guzman_neilan.py
@@ -109,15 +109,12 @@ def GuzmanNeilanH1div(ref_el, degree=2, reduced=False):
     """
     order = 0
     AS = AlfeldSorokina(ref_el, 2)
-    if reduced:
+    if reduced or ref_el.get_spatial_dimension() <= 2:
         order = 1
         # Only extract the div bubbles
         div_nodes = [i for i, node in enumerate(AS.dual_basis())
                      if len(node.deriv_dict) > 0]
         AS = RestrictedElement(AS, indices=div_nodes)
-    elif ref_el.get_spatial_dimension() <= 2:
-        # Quadratic bubbles are already included in 2D
-        return AS
     GN = GuzmanNeilanH1(ref_el, order=order)
     return NodalEnrichedElement(AS, GN)
 
diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py
new file mode 100644
index 000000000..5068a6b39
--- /dev/null
+++ b/FIAT/hu_zhang.py
@@ -0,0 +1,89 @@
+# -*- coding: utf-8 -*-
+"""Implementation of the Hu-Zhang finite elements."""
+
+# Copyright (C) 2024 by Francis Aznaran (University of Notre Dame)
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+
+from FIAT import finite_element, polynomial_set, dual_set
+from FIAT.check_format_variant import check_format_variant
+from FIAT.reference_element import TRIANGLE
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.functional import (ComponentPointEvaluation,
+                             PointwiseInnerProductEvaluation,
+                             TensorBidirectionalIntegralMoment,
+                             IntegralLegendreNormalNormalMoment,
+                             IntegralLegendreNormalTangentialMoment)
+
+
+class HuZhangDual(dual_set.DualSet):
+    def __init__(self, ref_el, degree, variant, qdegree):
+        top = ref_el.get_topology()
+        sd = ref_el.get_spatial_dimension()
+        shp = (sd, sd)
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+        nodes = []
+
+        # vertex dofs
+        for v in sorted(top[0]):
+            cur = len(nodes)
+            pt, = ref_el.make_points(0, v, degree)
+            nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt)
+                         for i in range(sd) for j in range(i, sd))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        # edge dofs
+        for entity in sorted(top[1]):
+            cur = len(nodes)
+            if variant == "point":
+                # nn and nt components evaluated at edge points
+                n = ref_el.compute_scaled_normal(entity)
+                t = ref_el.compute_edge_tangent(entity)
+                pts = ref_el.make_points(1, entity, degree)
+                nodes.extend(PointwiseInnerProductEvaluation(ref_el, n, s, pt)
+                             for pt in pts for s in (n, t))
+
+            elif variant == "integral":
+                # bidirectional nn and nt moments against P_{k-2}
+                moments = (IntegralLegendreNormalNormalMoment, IntegralLegendreNormalTangentialMoment)
+                nodes.extend(mu(ref_el, entity, order, qdegree + degree-2)
+                             for order in range(degree-1) for mu in moments)
+            entity_ids[1][entity].extend(range(cur, len(nodes)))
+
+        # interior dofs
+        cur = len(nodes)
+        if variant == "point":
+            # unique components evaluated at interior points
+            pts = ref_el.make_points(sd, 0, degree+1)
+            nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt)
+                         for pt in pts for i in range(sd) for j in range(i, sd))
+
+        elif variant == "integral":
+            # Moments of unique components against a basis for P_{k-2}
+            n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
+            Q = create_quadrature(ref_el, 2*degree-2)
+            P = polynomial_set.ONPolynomialSet(ref_el, degree-2, scale="L2 piola")
+            phis = P.tabulate(Q.get_points())[(0,)*sd]
+            nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
+                         for phi in phis for i in range(sd) for j in range(i, sd))
+
+        entity_ids[2][0].extend(range(cur, len(nodes)))
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class HuZhang(finite_element.CiarletElement):
+    """The definition of the Hu-Zhang element."""
+    def __init__(self, ref_el, degree=3, variant=None):
+        if degree < 3:
+            raise ValueError(f"{type(self).__name__} only defined for degree >= 3")
+        if ref_el.shape != TRIANGLE:
+            raise ValueError(f"{type(self).__name__} only defined on triangles")
+        variant, qdegree = check_format_variant(variant, degree)
+        poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree)
+        dual = HuZhangDual(ref_el, degree, variant, qdegree)
+        formdegree = ref_el.get_spatial_dimension() - 1
+        mapping = "double contravariant piola"
+        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py
index de4683bf4..a36907212 100644
--- a/FIAT/johnson_mercier.py
+++ b/FIAT/johnson_mercier.py
@@ -1,5 +1,5 @@
 from FIAT import finite_element, dual_set, macro, polynomial_set
-from FIAT.functional import IntegralMoment, FrobeniusIntegralMoment
+from FIAT.functional import TensorBidirectionalIntegralMoment
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.quadrature_schemes import create_quadrature
 import numpy
@@ -18,33 +18,30 @@ def __init__(self, ref_complex, degree, variant=None):
         nodes = []
 
         # Face dofs: bidirectional (nn and nt) Legendre moments
-        R = numpy.array([[0, 1], [-1, 0]])
         dim = sd - 1
+        R = numpy.array([[0, 1], [-1, 0]])
         ref_facet = ref_el.construct_subelement(dim)
         Qref = create_quadrature(ref_facet, 2*degree)
         P = polynomial_set.ONPolynomialSet(ref_facet, degree)
         phis = P.tabulate(Qref.get_points())[(0,) * dim]
-        for facet in sorted(top[dim]):
+        for f in sorted(top[dim]):
             cur = len(nodes)
-            Q = FacetQuadratureRule(ref_el, dim, facet, Qref)
-            thats = ref_el.compute_tangents(dim, facet)
+            Q = FacetQuadratureRule(ref_el, dim, f, Qref)
+            thats = ref_el.compute_tangents(dim, f)
             nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats)
             normal = nhat / Q.jacobian_determinant()
-
-            uvecs = (nhat, *thats)
-            comps = [numpy.outer(normal, uvec) for uvec in uvecs]
-            nodes.extend(FrobeniusIntegralMoment(ref_el, Q, comp[:, :, None] * phi[None, None, :])
-                         for phi in phis for comp in comps)
-            entity_ids[dim][facet].extend(range(cur, len(nodes)))
+            nodes.extend(TensorBidirectionalIntegralMoment(ref_el, normal, comp, Q, phi)
+                         for phi in phis for comp in (nhat, *thats))
+            entity_ids[dim][f].extend(range(cur, len(nodes)))
 
         cur = len(nodes)
-        if variant is None:
-            # Interior dofs: moments for each independent component
-            Q = create_quadrature(ref_complex, 2*degree-1)
-            P = polynomial_set.ONPolynomialSet(ref_el, degree-1)
-            phis = P.tabulate(Q.get_points())[(0,) * sd]
-            nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(i, j))
-                         for j in range(sd) for i in range(j+1) for phi in phis)
+        # Interior dofs: moments for each independent component
+        n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1])))
+        Q = create_quadrature(ref_complex, 2*degree-1)
+        P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola")
+        phis = P.tabulate(Q.get_points())[(0,) * sd]
+        nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi)
+                     for phi in phis for i in range(sd) for j in range(i, sd))
 
         entity_ids[sd][0].extend(range(cur, len(nodes)))
 
@@ -58,5 +55,6 @@ def __init__(self, ref_el, degree=1, variant=None):
         ref_complex = macro.AlfeldSplit(ref_el)
         poly_set = macro.HDivSymPolynomialSet(ref_complex, degree)
         dual = JohnsonMercierDualSet(ref_complex, degree, variant=variant)
+        formdegree = ref_el.get_spatial_dimension() - 1
         mapping = "double contravariant piola"
-        super().__init__(poly_set, dual, degree, mapping=mapping)
+        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py
index 0dbc76262..2c623c64f 100644
--- a/FIAT/mardal_tai_winther.py
+++ b/FIAT/mardal_tai_winther.py
@@ -8,154 +8,98 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 
-from FIAT.finite_element import CiarletElement
-from FIAT.dual_set import DualSet
-from FIAT.polynomial_set import ONPolynomialSet
+from FIAT import dual_set, expansions, finite_element, polynomial_set
 from FIAT.functional import (IntegralMomentOfNormalEvaluation,
                              IntegralMomentOfTangentialEvaluation,
                              IntegralLegendreNormalMoment,
                              IntegralMomentOfDivergence)
 
-from FIAT.quadrature import make_quadrature
+from FIAT.quadrature_schemes import create_quadrature
 
 
-def DivergenceDubinerMoments(cell, start_deg, stop_deg, comp_deg):
-    onp = ONPolynomialSet(cell, stop_deg)
-    Q = make_quadrature(cell, comp_deg)
+def DivergenceDubinerMoments(ref_el, start_deg, stop_deg, comp_deg):
+    sd = ref_el.get_spatial_dimension()
+    P = polynomial_set.ONPolynomialSet(ref_el, stop_deg)
+    Q = create_quadrature(ref_el, comp_deg + stop_deg)
 
-    pts = Q.get_points()
-    onp = onp.tabulate(pts, 0)[0, 0]
+    dim0 = expansions.polynomial_dimension(ref_el, start_deg-1)
+    dim1 = expansions.polynomial_dimension(ref_el, stop_deg)
+    indices = list(range(dim0, dim1))
+    phis = P.take(indices).tabulate(Q.get_points())[(0,)*sd]
+    for phi in phis:
+        yield IntegralMomentOfDivergence(ref_el, Q, phi)
 
-    ells = []
 
-    for ii in range((start_deg)*(start_deg+1)//2,
-                    (stop_deg+1)*(stop_deg+2)//2):
-        ells.append(IntegralMomentOfDivergence(cell, Q, onp[ii, :]))
+class MardalTaiWintherDual(dual_set.DualSet):
+    """Degrees of freedom for Mardal-Tai-Winther elements."""
+    def __init__(self, ref_el, degree):
+        sd = ref_el.get_spatial_dimension()
+        top = ref_el.get_topology()
 
-    return ells
+        if sd != 2:
+            raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2.")
 
+        if degree != 3:
+            raise ValueError("Mardal-Tai-Winther elements are only defined for degree 3.")
 
-class MardalTaiWintherDual(DualSet):
-    """Degrees of freedom for Mardal-Tai-Winther elements."""
-    def __init__(self, cell, degree):
-        dim = cell.get_spatial_dimension()
-        if not dim == 2:
-            raise ValueError("Mardal-Tai-Winther elements are only"
-                             "defined in dimension 2.")
-
-        if not degree == 3:
-            raise ValueError("Mardal-Tai-Winther elements are only defined"
-                             "for degree 3.")
-
-        # construct the degrees of freedoms
-        dofs = []               # list of functionals
-
-        # dof_ids[i][j] contains the indices of dofs that are associated with
-        # entity j in dim i
-        dof_ids = {}
-
-        # no vertex dof
-        dof_ids[0] = {i: [] for i in range(dim + 1)}
-
-        # edge dofs
-        (_dofs, _dof_ids) = self._generate_edge_dofs(cell, degree)
-        dofs.extend(_dofs)
-        dof_ids[1] = _dof_ids
-
-        # no cell dofs
-        dof_ids[2] = {}
-        dof_ids[2][0] = []
-
-        # extra dofs for enforcing div(v) constant over the cell and
-        # v.n linear on edges
-        (_dofs, _edge_dof_ids, _cell_dof_ids) = self._generate_constraint_dofs(cell, degree, len(dofs))
-        dofs.extend(_dofs)
-
-        for entity_id in range(3):
-            dof_ids[1][entity_id] = dof_ids[1][entity_id] + _edge_dof_ids[entity_id]
-
-        dof_ids[2][0] = dof_ids[2][0] + _cell_dof_ids
-
-        super(MardalTaiWintherDual, self).__init__(dofs, cell, dof_ids)
-
-    @staticmethod
-    def _generate_edge_dofs(cell, degree):
-        """Generate dofs on edges.
-        On each edge, let n be its normal.  We need to integrate
-        u.n and u.t against the first Legendre polynomial (constant)
-        and u.n against the second (linear).
-        """
-        dofs = []
-        dof_ids = {}
-        offset = 0
-        sd = 2
-
-        facet = cell.get_facet_element()
-        # Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1}
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
+
+        # no vertex dofs
+
+        # On each facet, let n be its normal.  We need to integrate
+        # u.n and u.t against the first Legendre polynomial (constant)
+        # and u.n against the second (linear).
+        facet = ref_el.get_facet_element()
+        # Facet nodes are \int_F v.n p ds where p \in P_{q-1}
         # degree is q - 1
-        Q = make_quadrature(facet, 6)
-        Pq = ONPolynomialSet(facet, 1)
-        Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(sd - 1))]
-        for f in range(3):
-            phi0 = Pq_at_qpts[0, :]
-            dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phi0, f))
-            dofs.append(IntegralMomentOfTangentialEvaluation(cell, Q, phi0, f))
-            phi1 = Pq_at_qpts[1, :]
-            dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phi1, f))
-
-            num_new_dofs = 3
-            dof_ids[f] = list(range(offset, offset + num_new_dofs))
-            offset += num_new_dofs
-
-        return (dofs, dof_ids)
-
-    @staticmethod
-    def _generate_constraint_dofs(cell, degree, offset):
-        """
-        Generate constraint dofs on the cell and edges
-        * div(v) must be constant on the cell.  Since v is a cubic and
-          div(v) is quadratic, we need the integral of div(v) against the
-          linear and quadratic Dubiner polynomials to vanish.
-          There are two linear and three quadratics, so these are five
-          constraints
-        * v.n must be linear on each edge.  Since v.n is cubic, we need
-          the integral of v.n against the cubic and quadratic Legendre
-          polynomial to vanish on each edge.
-
-        So we introduce functionals whose kernel describes this property,
-        as described in the FIAT paper.
-        """
-        dofs = []
-
-        edge_dof_ids = {}
-        for entity_id in range(3):
-            dofs += [IntegralLegendreNormalMoment(cell, entity_id, 2, 6),
-                     IntegralLegendreNormalMoment(cell, entity_id, 3, 6)]
-
-            edge_dof_ids[entity_id] = [offset, offset+1]
-            offset += 2
-
-        cell_dofs = DivergenceDubinerMoments(cell, 1, 2, 6)
-        dofs.extend(cell_dofs)
-        cell_dof_ids = list(range(offset, offset+len(cell_dofs)))
-
-        return (dofs, edge_dof_ids, cell_dof_ids)
-
-
-class MardalTaiWinther(CiarletElement):
+        Q = create_quadrature(facet, degree+1)
+        Pq = polynomial_set.ONPolynomialSet(facet, 1)
+        phis = Pq.tabulate(Q.get_points())[(0,)*(sd - 1)]
+        for f in sorted(top[sd-1]):
+            cur = len(nodes)
+            nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[0], f))
+            nodes.append(IntegralMomentOfTangentialEvaluation(ref_el, Q, phis[0], f))
+            nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[1], f))
+            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+
+        # Generate constraint nodes on the cell and facets
+        # * div(v) must be constant on the cell.  Since v is a cubic and
+        #   div(v) is quadratic, we need the integral of div(v) against the
+        #   linear and quadratic Dubiner polynomials to vanish.
+        #   There are two linear and three quadratics, so these are five
+        #   constraints
+        # * v.n must be linear on each facet.  Since v.n is cubic, we need
+        #   the integral of v.n against the cubic and quadratic Legendre
+        #   polynomial to vanish on each facet.
+
+        # So we introduce functionals whose kernel describes this property,
+        # as described in the FIAT paper.
+        start_order = 2
+        stop_order = 3
+        qdegree = degree + stop_order
+        for f in sorted(top[sd-1]):
+            cur = len(nodes)
+            nodes.extend(IntegralLegendreNormalMoment(ref_el, f, order, qdegree)
+                         for order in range(start_order, stop_order+1))
+            entity_ids[sd-1][f].extend(range(cur, len(nodes)))
+
+        cur = len(nodes)
+        nodes.extend(DivergenceDubinerMoments(ref_el, start_order-1, stop_order-1, degree))
+        entity_ids[sd][0].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class MardalTaiWinther(finite_element.CiarletElement):
     """The definition of the Mardal-Tai-Winther element.
     """
-    def __init__(self, cell, degree=3):
+    def __init__(self, ref_el, degree=3):
+        sd = ref_el.get_spatial_dimension()
         assert degree == 3, "Only defined for degree 3"
-        assert cell.get_spatial_dimension() == 2, "Only defined for dimension 2"
-        # polynomial space
-        Ps = ONPolynomialSet(cell, degree, (2,))
-
-        # degrees of freedom
-        Ls = MardalTaiWintherDual(cell, degree)
-
-        # mapping under affine transformation
+        assert sd == 2, "Only defined for dimension 2"
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd,))
+        dual = MardalTaiWintherDual(ref_el, degree)
+        formdegree = sd-1
         mapping = "contravariant piola"
-
-        super(MardalTaiWinther, self).__init__(Ps, Ls, degree,
-                                               mapping=mapping)
+        super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)
diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py
index 579ec3614..6210ef0b4 100644
--- a/FIAT/polynomial_set.py
+++ b/FIAT/polynomial_set.py
@@ -125,16 +125,13 @@ def __init__(self, ref_el, degree, shape=tuple(), **kwargs):
         if shape == tuple():
             coeffs = numpy.eye(num_members)
         else:
-            coeffs_shape = (num_members, *shape, num_exp_functions)
-            coeffs = numpy.zeros(coeffs_shape, "d")
-            # use functional's index_iterator function
-            cur_bf = 0
+            coeffs = numpy.zeros((num_members, *shape, num_exp_functions))
+            cur = 0
+            exp_bf = range(num_exp_functions)
             for idx in index_iterator(shape):
-                n = expansion_set.get_num_members(embedded_degree)
-                for exp_bf in range(n):
-                    cur_idx = (cur_bf, *idx, exp_bf)
-                    coeffs[cur_idx] = 1.0
-                    cur_bf += 1
+                cur_bf = range(cur, cur+num_exp_functions)
+                coeffs[(cur_bf, *idx, exp_bf)] = 1.0
+                cur += num_exp_functions
 
         super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs)
 
@@ -209,19 +206,49 @@ def __init__(self, ref_el, degree, size=None, **kwargs):
         embedded_degree = degree
 
         # set up coefficients for symmetric tensors
-        coeffs_shape = (num_members, *shape, num_exp_functions)
-        coeffs = numpy.zeros(coeffs_shape, "d")
-        cur_bf = 0
+        coeffs = numpy.zeros((num_members, *shape, num_exp_functions))
+        cur = 0
+        exp_bf = range(num_exp_functions)
         for i, j in index_iterator(shape):
+            if i > j:
+                continue
+            cur_bf = range(cur, cur+num_exp_functions)
+            coeffs[cur_bf, i, j, exp_bf] = 1.0
+            coeffs[cur_bf, j, i, exp_bf] = 1.0
+            cur += num_exp_functions
+
+        super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs)
+
+
+class TracelessTensorPolynomialSet(PolynomialSet):
+    """Constructs an orthonormal basis for traceless-tensor-valued
+    polynomials on a reference element.
+    """
+    def __init__(self, ref_el, degree, size=None, **kwargs):
+        expansion_set = expansions.ExpansionSet(ref_el, **kwargs)
+
+        sd = ref_el.get_spatial_dimension()
+        if size is None:
+            size = sd
+
+        shape = (size, size)
+        num_exp_functions = expansion_set.get_num_members(degree)
+        num_components = size * size - 1
+        num_members = num_components * num_exp_functions
+        embedded_degree = degree
+
+        # set up coefficients for traceless tensors
+        coeffs = numpy.zeros((num_members, *shape, num_exp_functions))
+        cur = 0
+        exp_bf = range(num_exp_functions)
+        for i, j in index_iterator(shape):
+            if i == size-1 and j == size-1:
+                continue
+            cur_bf = range(cur, cur+num_exp_functions)
+            coeffs[cur_bf, i, j, exp_bf] = 1.0
             if i == j:
-                for exp_bf in range(num_exp_functions):
-                    coeffs[cur_bf, i, j, exp_bf] = 1.0
-                    cur_bf += 1
-            elif i < j:
-                for exp_bf in range(num_exp_functions):
-                    coeffs[cur_bf, i, j, exp_bf] = 1.0
-                    coeffs[cur_bf, j, i, exp_bf] = 1.0
-                    cur_bf += 1
+                coeffs[cur_bf, -1, -1, exp_bf] = -1.0
+            cur += num_exp_functions
 
         super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs)
 
diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index f09751578..742df82e3 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -570,7 +570,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
 
     def compute_bubble(self, points, entity=None):
         """Returns the lowest-order bubble on an entity evaluated at the given
-        points on the entity."""
+        points on the cell."""
         return numpy.prod(self.compute_barycentric_coordinates(points, entity), axis=1)
 
     def distance_to_point_l1(self, points, entity=None, rescale=False):
diff --git a/test/unit/test_awc.py b/test/unit/test_awc.py
index f40690c48..b96bc86be 100644
--- a/test/unit/test_awc.py
+++ b/test/unit/test_awc.py
@@ -1,11 +1,11 @@
 import numpy as np
-from FIAT import ufc_simplex, ArnoldWinther, make_quadrature, expansions
+from FIAT import ufc_simplex, ArnoldWinther, create_quadrature, expansions
 
 
 def test_dofs():
     line = ufc_simplex(1)
     T = ufc_simplex(2)
-    T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)])
+    T.vertices = ((0.0, 0.0), (1.0, 0.25), (-0.75, 1.1))
     AW = ArnoldWinther(T, 3)
 
     # check Kronecker property at vertices
@@ -20,7 +20,7 @@ def test_dofs():
                 assert np.allclose(vert_vals[3*i+j, :, :, (i+k) % 3], np.zeros((2, 2)))
 
     # check edge moments
-    Qline = make_quadrature(line, 6)
+    Qline = create_quadrature(line, 6)
 
     linebfs = expansions.LineExpansionSet(line)
     linevals = linebfs.tabulate(1, Qline.pts)
@@ -73,15 +73,21 @@ def test_dofs():
                 assert np.allclose(ntmoments[bf, :], np.zeros(2))
 
     # check internal dofs
-    Q = make_quadrature(T, 6)
+    ns = list(map(T.compute_scaled_normal, range(3)))
+    Q = create_quadrature(T, 3)
     qpvals = AW.tabulate(0, Q.pts)[(0, 0)]
-    const_moms = qpvals @ Q.wts
-    assert np.allclose(const_moms[:21], np.zeros((21, 2, 2)))
-    assert np.allclose(const_moms[24:], np.zeros((6, 2, 2)))
-    assert np.allclose(const_moms[21:24, 0, 0], np.asarray([1, 0, 0]))
-    assert np.allclose(const_moms[21:24, 0, 1], np.asarray([0, 1, 0]))
-    assert np.allclose(const_moms[21:24, 1, 0], np.asarray([0, 1, 0]))
-    assert np.allclose(const_moms[21:24, 1, 1], np.asarray([0, 0, 1]))
+    const_moms = qpvals @ Q.wts / T.volume()
+    nn_moms = const_moms.copy()
+    for j in range(2):
+        for i in range(2):
+            comp = np.outer(ns[i+1], ns[j+1])
+            nn_moms[:, i, j] = np.tensordot(const_moms, comp, ((1, 2), (0, 1)))
+    assert np.allclose(nn_moms[:21], np.zeros((21, 2, 2)))
+    assert np.allclose(nn_moms[24:], np.zeros((6, 2, 2)))
+    assert np.allclose(nn_moms[21:24, 0, 0], np.asarray([1, 0, 0]))
+    assert np.allclose(nn_moms[21:24, 0, 1], np.asarray([0, 1, 0]))
+    assert np.allclose(nn_moms[21:24, 1, 0], np.asarray([0, 1, 0]))
+    assert np.allclose(nn_moms[21:24, 1, 1], np.asarray([0, 0, 1]))
 
 
 def frob(a, b):
@@ -90,11 +96,11 @@ def frob(a, b):
 
 def test_projection():
     T = ufc_simplex(2)
-    T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.0), (0.5, 2.1)])
+    T.vertices = ((0.0, 0.0), (1.0, 0.0), (0.5, 2.1))
 
     AW = ArnoldWinther(T, 3)
 
-    Q = make_quadrature(T, 4)
+    Q = create_quadrature(T, 6)
     qpts = np.asarray(Q.pts)
     qwts = np.asarray(Q.wts)
     nqp = len(Q.wts)
diff --git a/test/unit/test_awnc.py b/test/unit/test_awnc.py
index f7a28cac8..c770791ee 100644
--- a/test/unit/test_awnc.py
+++ b/test/unit/test_awnc.py
@@ -5,7 +5,7 @@
 def test_dofs():
     line = ufc_simplex(1)
     T = ufc_simplex(2)
-    T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)])
+    T.vertices = ((0.0, 0.0), (1.0, 0.25), (-0.75, 1.1))
     AW = ArnoldWintherNC(T, 2)
 
     Qline = make_quadrature(line, 6)
@@ -61,12 +61,18 @@ def test_dofs():
                 assert np.allclose(ntmoments[bf, :], np.zeros(2), atol=1.e-7)
 
     # check internal dofs
+    ns = list(map(T.compute_scaled_normal, range(3)))
     Q = make_quadrature(T, 6)
     qpvals = AW.tabulate(0, Q.pts)[(0, 0)]
-    const_moms = qpvals @ Q.wts
-    assert np.allclose(const_moms[:12], np.zeros((12, 2, 2)))
-    assert np.allclose(const_moms[15:], np.zeros((3, 2, 2)))
-    assert np.allclose(const_moms[12:15, 0, 0], np.asarray([1, 0, 0]))
-    assert np.allclose(const_moms[12:15, 0, 1], np.asarray([0, 1, 0]))
-    assert np.allclose(const_moms[12:15, 1, 0], np.asarray([0, 1, 0]))
-    assert np.allclose(const_moms[12:15, 1, 1], np.asarray([0, 0, 1]))
+    const_moms = qpvals @ Q.wts / T.volume()
+    nn_moms = const_moms.copy()
+    for j in range(2):
+        for i in range(2):
+            comp = np.outer(ns[i+1], ns[j+1])
+            nn_moms[:, i, j] = np.tensordot(const_moms, comp, ((1, 2), (0, 1)))
+    assert np.allclose(nn_moms[:12], np.zeros((12, 2, 2)))
+    assert np.allclose(nn_moms[15:], np.zeros((3, 2, 2)))
+    assert np.allclose(nn_moms[12:15, 0, 0], np.asarray([1, 0, 0]))
+    assert np.allclose(nn_moms[12:15, 0, 1], np.asarray([0, 1, 0]))
+    assert np.allclose(nn_moms[12:15, 1, 0], np.asarray([0, 1, 0]))
+    assert np.allclose(nn_moms[12:15, 1, 1], np.asarray([0, 0, 1]))
diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py
index f074d61b0..806e74911 100644
--- a/test/unit/test_fiat.py
+++ b/test/unit/test_fiat.py
@@ -35,6 +35,8 @@
 from FIAT.regge import Regge                                    # noqa: F401
 from FIAT.hdiv_trace import HDivTrace, map_to_reference_facet   # noqa: F401
 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson  # noqa: F401
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind  # noqa: F401
+from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind  # noqa: F401
 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini  # noqa: F401
 from FIAT.gauss_legendre import GaussLegendre                   # noqa: F401
 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre    # noqa: F401
@@ -42,6 +44,9 @@
 from FIAT.tensor_product import TensorProductElement            # noqa: F401
 from FIAT.tensor_product import FlattenedDimensions             # noqa: F401
 from FIAT.hdivcurl import Hdiv, Hcurl                           # noqa: F401
+from FIAT.mardal_tai_winther import MardalTaiWinther            # noqa: F401
+from FIAT.arnold_winther import ArnoldWinther, ArnoldWintherNC  # noqa: F401
+from FIAT.hu_zhang import HuZhang                               # noqa: F401
 from FIAT.bernardi_raugel import BernardiRaugel                 # noqa: F401
 from FIAT.argyris import Argyris                                # noqa: F401
 from FIAT.hermite import CubicHermite                           # noqa: F401
@@ -270,6 +275,18 @@ def __init__(self, a, b):
     "HellanHerrmannJohnson(S, 2)",
     "HellanHerrmannJohnson(T, 1, variant='point')",
     "HellanHerrmannJohnson(S, 1, variant='point')",
+    "GopalakrishnanLedererSchoberlFirstKind(T, 1)",
+    "GopalakrishnanLedererSchoberlFirstKind(T, 2)",
+    "GopalakrishnanLedererSchoberlFirstKind(T, 3)",
+    "GopalakrishnanLedererSchoberlFirstKind(S, 1)",
+    "GopalakrishnanLedererSchoberlFirstKind(S, 2)",
+    "GopalakrishnanLedererSchoberlFirstKind(S, 3)",
+    "GopalakrishnanLedererSchoberlSecondKind(T, 0)",
+    "GopalakrishnanLedererSchoberlSecondKind(T, 1)",
+    "GopalakrishnanLedererSchoberlSecondKind(T, 2)",
+    "GopalakrishnanLedererSchoberlSecondKind(S, 0)",
+    "GopalakrishnanLedererSchoberlSecondKind(S, 1)",
+    "GopalakrishnanLedererSchoberlSecondKind(S, 2)",
     "BrezziDouglasFortinMarini(T, 2)",
     "GaussLegendre(I, 0)",
     "GaussLegendre(I, 1)",
@@ -315,6 +332,13 @@ def __init__(self, a, b):
     "Morley(T)",
     "BernardiRaugel(T)",
     "BernardiRaugel(S)",
+    "MardalTaiWinther(T, 3)",
+    "ArnoldWintherNC(T, 2)",
+    "ArnoldWinther(T, 3)",
+    "HuZhang(T, 3)",
+    "HuZhang(T, 4)",
+    "HuZhang(T, 3, 'point')",
+    "HuZhang(T, 4, 'point')",
 
     # Macroelements
     "Lagrange(T, 1, 'iso')",
diff --git a/test/unit/test_gopalakrishnan_lederer_schoberl.py b/test/unit/test_gopalakrishnan_lederer_schoberl.py
new file mode 100644
index 000000000..f78c8da34
--- /dev/null
+++ b/test/unit/test_gopalakrishnan_lederer_schoberl.py
@@ -0,0 +1,77 @@
+import pytest
+import numpy
+
+from FIAT import (GopalakrishnanLedererSchoberlFirstKind,
+                  GopalakrishnanLedererSchoberlSecondKind)
+from FIAT.reference_element import ufc_simplex
+from FIAT.expansions import polynomial_dimension
+from FIAT.polynomial_set import ONPolynomialSet
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.quadrature import FacetQuadratureRule
+
+
+@pytest.fixture(params=("T", "S"))
+def cell(request):
+    dim = {"I": 1, "T": 2, "S": 3}[request.param]
+    return ufc_simplex(dim)
+
+
+@pytest.mark.parametrize("degree", (1, 2, 3))
+@pytest.mark.parametrize("kind", (1, 2))
+def test_gls_bubbles(kind, cell, degree):
+    if kind == 1:
+        element = GopalakrishnanLedererSchoberlFirstKind
+    else:
+        element = GopalakrishnanLedererSchoberlSecondKind
+    fe = element(cell, degree)
+    sd = cell.get_spatial_dimension()
+    facet_el = cell.construct_subelement(sd-1)
+    poly_set = fe.get_nodal_basis()
+
+    # test dimension of constrained space
+    dimPkm1 = polynomial_dimension(facet_el, degree-1)
+    dimPkp1 = polynomial_dimension(facet_el, degree+1)
+    dimPk = polynomial_dimension(facet_el, degree)
+    if kind == 1:
+        constraints = dimPk - dimPkm1
+    else:
+        constraints = 0
+    expected = (sd**2-1)*(polynomial_dimension(cell, degree) - constraints)
+    assert poly_set.get_num_members() == expected
+
+    # test dimension of the bubbles
+    entity_dofs = fe.entity_dofs()
+    bubbles = poly_set.take(entity_dofs[sd][0])
+    expected = (sd**2-1)*polynomial_dimension(cell, degree-1)
+    assert bubbles.get_num_members() == expected
+
+    top = cell.get_topology()
+    Qref = create_quadrature(facet_el, 2*degree+1)
+    Pk = ONPolynomialSet(facet_el, degree+1)
+    if kind == 1:
+        start, stop = dimPkm1, dimPkp1
+    else:
+        start, stop = dimPk, dimPkp1
+    PkH = Pk.take(list(range(start, stop)))
+    PkH_at_qpts = PkH.tabulate(Qref.get_points())[(0,)*(sd-1)]
+    weights = numpy.transpose(numpy.multiply(PkH_at_qpts, Qref.get_weights()))
+    for facet in top[sd-1]:
+        n = cell.compute_scaled_normal(facet)
+        rts = cell.compute_tangents(sd-1, facet)
+        Q = FacetQuadratureRule(cell, sd-1, facet, Qref)
+        qpts, qwts = Q.get_points(), Q.get_weights()
+
+        # test the degree of normal-tangential components
+        phi_at_pts = fe.tabulate(0, qpts)[(0,) * sd]
+        for t in rts:
+            nt = numpy.outer(t, n)
+            phi_nt = numpy.tensordot(nt, phi_at_pts, axes=((0, 1), (1, 2)))
+            assert numpy.allclose(numpy.dot(phi_nt, weights), 0)
+
+        # test the support of the normal-tangential bubble
+        phi_at_pts = bubbles.tabulate(qpts)[(0,) * sd]
+        for t in rts:
+            nt = numpy.outer(t, n)
+            phi_nt = numpy.tensordot(nt, phi_at_pts, axes=((0, 1), (1, 2)))
+            norms = numpy.dot(phi_nt**2, qwts)
+            assert numpy.allclose(norms, 0)