From 90978c9d386cb0711250069f6225055f1db0d8cd Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Fri, 27 Jun 2025 11:09:23 +0100
Subject: [PATCH 01/11] Bubble Lagrange variants

---
 FIAT/bubble.py           | 12 ++++++------
 finat/element_factory.py |  2 +-
 finat/fiat_elements.py   |  8 ++++----
 3 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/FIAT/bubble.py b/FIAT/bubble.py
index 0de7b2566..6ecfeb826 100644
--- a/FIAT/bubble.py
+++ b/FIAT/bubble.py
@@ -14,8 +14,8 @@
 class CodimBubble(RestrictedElement):
     """Bubbles of a certain codimension."""
 
-    def __init__(self, ref_el, degree, codim):
-        element = Lagrange(ref_el, degree)
+    def __init__(self, ref_el, degree, codim, variant=None):
+        element = Lagrange(ref_el, degree, variant=variant)
 
         cell_dim = ref_el.get_dimension()
         assert cell_dim == max(element.entity_dofs().keys())
@@ -29,12 +29,12 @@ def __init__(self, ref_el, degree, codim):
 class Bubble(CodimBubble):
     """The bubble finite element: the dofs of the Lagrange FE in the interior of the cell"""
 
-    def __init__(self, ref_el, degree):
-        super().__init__(ref_el, degree, codim=0)
+    def __init__(self, ref_el, degree, variant=None):
+        super().__init__(ref_el, degree, codim=0, variant=variant)
 
 
 class FacetBubble(CodimBubble):
     """The facet bubble finite element: the dofs of the Lagrange FE in the interior of the facets"""
 
-    def __init__(self, ref_el, degree):
-        super().__init__(ref_el, degree, codim=1)
+    def __init__(self, ref_el, degree, variant=None):
+        super().__init__(ref_el, degree, codim=1, variant=variant)
diff --git a/finat/element_factory.py b/finat/element_factory.py
index 8e038dccd..51c4d3d66 100644
--- a/finat/element_factory.py
+++ b/finat/element_factory.py
@@ -219,7 +219,7 @@ def convert_finiteelement(element, **kwargs):
             lmbda = finat.DiscontinuousLagrange
             finat_kwargs["variant"] = kind
 
-    elif element.family() == "HDiv Trace":
+    elif element.family() in {"HDiv Trace", "Bubble", "FacetBubble"}:
         finat_kwargs["variant"] = kind
 
     elif element.variant() is not None:
diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py
index c4799b1ae..3677ef965 100644
--- a/finat/fiat_elements.py
+++ b/finat/fiat_elements.py
@@ -348,13 +348,13 @@ def __init__(self, cell, degree):
 
 
 class Bubble(ScalarFiatElement):
-    def __init__(self, cell, degree):
-        super().__init__(FIAT.Bubble(cell, degree))
+    def __init__(self, cell, degree, variant=None):
+        super().__init__(FIAT.Bubble(cell, degree, variant=variant))
 
 
 class FacetBubble(ScalarFiatElement):
-    def __init__(self, cell, degree):
-        super().__init__(FIAT.FacetBubble(cell, degree))
+    def __init__(self, cell, degree, variant=None):
+        super().__init__(FIAT.FacetBubble(cell, degree, variant=variant))
 
 
 class CrouzeixRaviart(ScalarFiatElement):

From 0499c1b7fb6578a05dbe84f57a022cad5c07ed46 Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Fri, 27 Jun 2025 17:10:01 +0100
Subject: [PATCH 02/11] HDivTrace: tabulate on multiple facets

---
 FIAT/bubble.py     |  4 ++-
 FIAT/hdiv_trace.py | 82 +++++++++++++++++++++-------------------------
 2 files changed, 41 insertions(+), 45 deletions(-)

diff --git a/FIAT/bubble.py b/FIAT/bubble.py
index 6ecfeb826..8d685b46e 100644
--- a/FIAT/bubble.py
+++ b/FIAT/bubble.py
@@ -7,6 +7,7 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 from FIAT.lagrange import Lagrange
+from FIAT.hierarchical import IntegratedLegendre
 from FIAT.restricted import RestrictedElement
 from itertools import chain
 
@@ -15,7 +16,8 @@ class CodimBubble(RestrictedElement):
     """Bubbles of a certain codimension."""
 
     def __init__(self, ref_el, degree, codim, variant=None):
-        element = Lagrange(ref_el, degree, variant=variant)
+        CG = IntegratedLegendre if variant == "integral" else Lagrange
+        element = CG(ref_el, degree, variant=variant)
 
         cell_dim = ref_el.get_dimension()
         assert cell_dim == max(element.entity_dofs().keys())
diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py
index 264205a3e..dab911625 100644
--- a/FIAT/hdiv_trace.py
+++ b/FIAT/hdiv_trace.py
@@ -5,8 +5,10 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 import numpy as np
+from collections import defaultdict
 from FIAT.barycentric_interpolation import get_lagrange_points
 from FIAT.discontinuous_lagrange import DiscontinuousLagrange
+from FIAT.hierarchical import Legendre
 from FIAT.dual_set import DualSet
 from FIAT.finite_element import FiniteElement
 from FIAT.functional import PointEvaluation
@@ -99,6 +101,7 @@ def __init__(self, ref_el, degree, variant=None):
             nf = element.space_dimension()
             num_facets = len(topology[facet_dim])
 
+            # FIXME non PointEvaluation dofs
             facet_pts = get_lagrange_points(element.dual_basis())
             for i in range(num_facets):
                 entity_dofs[facet_dim][i] = list(range(offset, offset + nf))
@@ -161,6 +164,7 @@ def tabulate(self, order, points, entity=None):
         """
         sd = self.ref_el.get_spatial_dimension()
         facet_sd = sd - 1
+        evalkey = (0,) * sd
 
         # Initializing dictionary with zeros
         phivals = {}
@@ -168,12 +172,13 @@ def tabulate(self, order, points, entity=None):
             alphas = mis(sd, i)
             for alpha in alphas:
                 phivals[alpha] = np.zeros(shape=(self.space_dimension(), len(points)))
-
-        evalkey = (0,) * sd
+                if alpha != evalkey:
+                    # If asking for gradient evaluations, insert TraceError in gradient slots
+                    phivals[alpha] = TraceError("Gradients on trace elements are not well-defined.")
 
         # If entity is None, identify facet using numerical tolerance and
         # return the tabulated values
-        if entity is None:
+        if entity is None or entity == (sd, 0):
             # NOTE: Numerical approximation of the facet id is currently only
             # implemented for simplex reference cells.
             if self.ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON]:
@@ -185,25 +190,26 @@ def tabulate(self, order, points, entity=None):
             # Attempt to identify which facet (if any) the given points are on
             vertices = self.ref_el.vertices
             coordinates = barycentric_coordinates(points, vertices)
-            unique_facet, success = extract_unique_facet(coordinates)
+            bins, success = extract_facets(coordinates)
 
             # If not successful, return NaNs
             if not success:
                 for key in phivals:
-                    phivals[key] = np.full(shape=(self.space_dimension(), len(points)), fill_value=np.nan)
-
-                return phivals
+                    phivals[key] = np.nan(shape=(self.space_dimension(), len(points)))
 
             # Otherwise, extract non-zero values and insertion indices
-            else:
+            for facet, ids in bins.items():
                 # Map points to the reference facet
-                new_points = map_to_reference_facet(points, vertices, unique_facet)
+                new_points = map_to_reference_facet(points[ids], vertices, facet)
 
                 # Retrieve values by tabulating the DG element
                 element = self.dg_elements[facet_sd]
                 nf = element.space_dimension()
-                nonzerovals, = element.tabulate(order, new_points).values()
-                indices = slice(nf * unique_facet, nf * (unique_facet + 1))
+                nonzerovals = element.tabulate(order, new_points)[(0,)*facet_sd]
+                indices = slice(nf * facet, nf * (facet + 1))
+
+                # Insert non-zero values in appropriate place
+                phivals[evalkey][indices, ids] = nonzerovals
 
         else:
             entity_dim, _ = entity
@@ -215,8 +221,6 @@ def tabulate(self, order, points, entity=None):
                     msg = "The HDivTrace element can only be tabulated on facets."
                     phivals[key] = TraceError(msg)
 
-                return phivals
-
             else:
                 # Retrieve function evaluations (order = 0 case)
                 offset = 0
@@ -235,16 +239,8 @@ def tabulate(self, order, points, entity=None):
 
                         offset += nf
 
-        # If asking for gradient evaluations, insert TraceError in
-        # gradient slots
-        if order > 0:
-            msg = "Gradients on trace elements are not well-defined."
-            for key in phivals:
-                if key != evalkey:
-                    phivals[key] = TraceError(msg)
-
-        # Insert non-zero values in appropriate place
-        phivals[evalkey][indices, :] = nonzerovals
+                # Insert non-zero values in appropriate place
+                phivals[evalkey][indices] = nonzerovals
 
         return phivals
 
@@ -270,13 +266,14 @@ def construct_dg_element(ref_el, degree, variant):
     """Constructs a discontinuous galerkin element of a given degree
     on a particular reference cell.
     """
+    DG = Legendre if variant == "integral" else DiscontinuousLagrange
     if ref_el.get_shape() in [LINE, TRIANGLE]:
-        dg_element = DiscontinuousLagrange(ref_el, degree, variant)
+        dg_element = DG(ref_el, degree, variant)
 
     # Quadrilateral facets could be on a FiredrakeQuadrilateral.
     # In this case, we treat this as an interval x interval cell:
     elif ref_el.get_shape() == QUADRILATERAL:
-        dg_line = DiscontinuousLagrange(ufc_simplex(1), degree, variant)
+        dg_line = DG(ufc_simplex(1), degree, variant)
         dg_element = TensorProductElement(dg_line, dg_line)
 
     # This handles the more general case for facets:
@@ -303,29 +300,26 @@ def construct_dg_element(ref_el, degree, variant):
 
 
 # The following functions are credited to Marie E. Rognes:
-def extract_unique_facet(coordinates, tolerance=epsilon):
-    """Determines whether a set of points (described in its barycentric coordinates)
-    are all on one of the facet sub-entities, and return the particular facet and
-    whether the search has been successful.
+def extract_facets(coordinates, tolerance=epsilon):
+    """Determines whether a set of points (described in barycentric coordinates)
+    are on facet sub-entities, and return a dict mapping facets to
+    point indices and whether the search has been successful.
 
     :arg coordinates: A set of points described in barycentric coordinates.
     :arg tolerance: A fixed tolerance for geometric identifications.
     """
-    facets = []
-    for c in coordinates:
-        on_facet = set([i for (i, l) in enumerate(c) if abs(l) < tolerance])
-        facets += [on_facet]
-
-    unique_facet = facets[0]
-    for f in facets:
-        unique_facet = unique_facet & f
-
-    # Handle coordinates not on facets
-    if len(unique_facet) != 1:
-        return (None, False)
-
-    # If we have a unique facet, return it and success
-    return (unique_facet.pop(), True)
+    bins = defaultdict(list)
+    for ipt, c in enumerate(coordinates):
+        on_facet = set(i for (i, l) in enumerate(c) if abs(l) < tolerance)
+        try:
+            f, = on_facet
+        except ValueError:
+            # Handle coordinates not on facets
+            return ({}, False)
+        bins[f].append(ipt)
+
+    # If all points are on facets, return bins and success
+    return (bins, True)
 
 
 def barycentric_coordinates(points, vertices):

From 3ae45fb5c18d14d824a519c28bb089dbc8b9fc5f Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Fri, 27 Jun 2025 21:22:14 +0100
Subject: [PATCH 03/11] HDiv Trace integral variant

---
 FIAT/hdiv_trace.py | 28 +++++++++++++++++++---------
 1 file changed, 19 insertions(+), 9 deletions(-)

diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py
index dab911625..804ad294b 100644
--- a/FIAT/hdiv_trace.py
+++ b/FIAT/hdiv_trace.py
@@ -13,6 +13,7 @@
 from FIAT.finite_element import FiniteElement
 from FIAT.functional import PointEvaluation
 from FIAT.polynomial_set import mis
+from FIAT.quadrature import FacetQuadratureRule
 from FIAT.reference_element import (ufc_simplex, POINT,
                                     LINE, QUADRILATERAL,
                                     TRIANGLE, TETRAHEDRON,
@@ -95,28 +96,22 @@ def __init__(self, ref_el, degree, variant=None):
         # Compute the dof numbering for all facet entities
         # and extract nodes
         offset = 0
-        pts = []
+        nodes = []
         for facet_dim in sorted(dg_elements):
             element = dg_elements[facet_dim]
             nf = element.space_dimension()
             num_facets = len(topology[facet_dim])
 
-            # FIXME non PointEvaluation dofs
-            facet_pts = get_lagrange_points(element.dual_basis())
             for i in range(num_facets):
                 entity_dofs[facet_dim][i] = list(range(offset, offset + nf))
                 offset += nf
-
-                # Run over nodes and collect the points for point evaluations
-                transform = ref_el.get_entity_transform(facet_dim, i)
-                pts.extend(transform(facet_pts))
+                nodes.extend(transform_nodes(element.dual_basis(), ref_el, facet_dim, i))
 
         # Setting up dual basis - only point evaluations
-        nodes = [PointEvaluation(ref_el, pt) for pt in pts]
         dual = DualSet(nodes, ref_el, entity_dofs)
 
         # Degree of the element
-        deg = max([e.degree() for e in dg_elements.values()])
+        deg = max(e.degree() for e in dg_elements.values())
 
         super().__init__(ref_el, dual, order=deg,
                          formdegree=facet_sd,
@@ -299,6 +294,21 @@ def construct_dg_element(ref_el, degree, variant):
     return dg_element
 
 
+def transform_nodes(ells, ref_el, facet_dim, facet_id):
+    """Map functionals into a given facet."""
+    try:
+        facet_pts = get_lagrange_points(ells)
+        transform = ref_el.get_entity_transform(facet_dim, facet_id)
+        pts = transform(facet_pts)
+        for pt in pts:
+            yield PointEvaluation(ref_el, pt)
+    except ValueError:
+        Q_ref, = set(ell.Q for ell in ells)
+        Q = FacetQuadratureRule(ref_el, facet_dim, facet_id, Q_ref)
+        for ell in ells:
+            yield type(ell)(ref_el, Q, ell.f_at_qpts)
+
+
 # The following functions are credited to Marie E. Rognes:
 def extract_facets(coordinates, tolerance=epsilon):
     """Determines whether a set of points (described in barycentric coordinates)

From eb2bd8e035ad6dac82d5ab394daefae759ae6865 Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Fri, 27 Jun 2025 23:47:36 +0100
Subject: [PATCH 04/11] Tests

---
 FIAT/hdiv_trace.py               | 56 +++++++++++++++-----------------
 test/FIAT/unit/test_fiat.py      |  9 ++++-
 test/FIAT/unit/test_hdivtrace.py |  5 +--
 3 files changed, 38 insertions(+), 32 deletions(-)

diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py
index 804ad294b..54b008844 100644
--- a/FIAT/hdiv_trace.py
+++ b/FIAT/hdiv_trace.py
@@ -11,7 +11,7 @@
 from FIAT.hierarchical import Legendre
 from FIAT.dual_set import DualSet
 from FIAT.finite_element import FiniteElement
-from FIAT.functional import PointEvaluation
+from FIAT.functional import IntegralMoment, PointEvaluation
 from FIAT.polynomial_set import mis
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.reference_element import (ufc_simplex, POINT,
@@ -76,22 +76,18 @@ def __init__(self, ref_el, degree, variant=None):
             if isinstance(degree, tuple):
                 raise ValueError("Must have a tensor product cell if providing multiple degrees")
 
-        # Initialize entity dofs and construct the DG elements
-        # for the facets
+        # Initialize entity dofs
         facet_sd = sd - 1
-        dg_elements = {}
-        entity_dofs = {}
         topology = ref_el.get_topology()
-        for top_dim, entities in topology.items():
-            cell = ref_el.construct_subelement(top_dim)
-            entity_dofs[top_dim] = {}
+        entity_dofs = {dim: {entity: [] for entity in topology[dim]} for dim in topology}
 
-            # We have a facet entity!
-            if cell.get_spatial_dimension() == facet_sd:
-                dg_elements[top_dim] = construct_dg_element(cell, degree, variant)
-            # Initialize
-            for entity in entities:
-                entity_dofs[top_dim][entity] = []
+        as_int = lambda x: sum(x) if isinstance(x, tuple) else x
+        # construct the DG element for the facets
+        dg_elements = {}
+        for dim in topology:
+            if as_int(dim) == facet_sd:
+                cell = ref_el.construct_subelement(dim)
+                dg_elements[dim] = construct_dg_element(cell, degree, variant)
 
         # Compute the dof numbering for all facet entities
         # and extract nodes
@@ -185,17 +181,21 @@ def tabulate(self, order, points, entity=None):
             # Attempt to identify which facet (if any) the given points are on
             vertices = self.ref_el.vertices
             coordinates = barycentric_coordinates(points, vertices)
-            bins, success = extract_facets(coordinates)
+            facet_to_pts, success = extract_facets(coordinates)
 
             # If not successful, return NaNs
             if not success:
                 for key in phivals:
-                    phivals[key] = np.nan(shape=(self.space_dimension(), len(points)))
+                    if entity is None:
+                        phivals[key].fill(np.nan)
+                    else:
+                        msg = "The HDivTrace element can only be tabulated on facets."
+                        phivals[key] = TraceError(msg)
 
             # Otherwise, extract non-zero values and insertion indices
-            for facet, ids in bins.items():
+            for facet, ipts in facet_to_pts.items():
                 # Map points to the reference facet
-                new_points = map_to_reference_facet(points[ids], vertices, facet)
+                new_points = map_to_reference_facet(points[ipts], vertices, facet)
 
                 # Retrieve values by tabulating the DG element
                 element = self.dg_elements[facet_sd]
@@ -204,7 +204,7 @@ def tabulate(self, order, points, entity=None):
                 indices = slice(nf * facet, nf * (facet + 1))
 
                 # Insert non-zero values in appropriate place
-                phivals[evalkey][indices, ids] = nonzerovals
+                phivals[evalkey][indices, ipts] = nonzerovals
 
         else:
             entity_dim, _ = entity
@@ -222,16 +222,14 @@ def tabulate(self, order, points, entity=None):
                 for facet_dim in sorted(self.dg_elements):
                     element = self.dg_elements[facet_dim]
                     nf = element.space_dimension()
-                    num_facets = len(self.ref_el.get_topology()[facet_dim])
 
                     # Loop over the number of facets until we find a facet
                     # with matching dimension and id
-                    for i in range(num_facets):
+                    for i in self.ref_el.get_topology()[facet_dim]:
                         # Found it! Grab insertion indices
                         if (facet_dim, i) == entity:
-                            nonzerovals, = element.tabulate(0, points).values()
+                            nonzerovals = element.tabulate(0, points)[(0,)*facet_sd]
                             indices = slice(offset, offset + nf)
-
                         offset += nf
 
                 # Insert non-zero values in appropriate place
@@ -306,19 +304,19 @@ def transform_nodes(ells, ref_el, facet_dim, facet_id):
         Q_ref, = set(ell.Q for ell in ells)
         Q = FacetQuadratureRule(ref_el, facet_dim, facet_id, Q_ref)
         for ell in ells:
-            yield type(ell)(ref_el, Q, ell.f_at_qpts)
+            yield IntegralMoment(ref_el, Q, ell.f_at_qpts)
 
 
 # The following functions are credited to Marie E. Rognes:
 def extract_facets(coordinates, tolerance=epsilon):
     """Determines whether a set of points (described in barycentric coordinates)
-    are on facet sub-entities, and return a dict mapping facets to
+    are all on facet sub-entities, and return a dict mapping facets to
     point indices and whether the search has been successful.
 
     :arg coordinates: A set of points described in barycentric coordinates.
     :arg tolerance: A fixed tolerance for geometric identifications.
     """
-    bins = defaultdict(list)
+    facet_to_pts = defaultdict(list)
     for ipt, c in enumerate(coordinates):
         on_facet = set(i for (i, l) in enumerate(c) if abs(l) < tolerance)
         try:
@@ -326,10 +324,10 @@ def extract_facets(coordinates, tolerance=epsilon):
         except ValueError:
             # Handle coordinates not on facets
             return ({}, False)
-        bins[f].append(ipt)
+        facet_to_pts[f].append(ipt)
 
-    # If all points are on facets, return bins and success
-    return (bins, True)
+    # If all points are on facets, return indices and success
+    return (facet_to_pts, True)
 
 
 def barycentric_coordinates(points, vertices):
diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py
index ffb192afc..745518a75 100644
--- a/test/FIAT/unit/test_fiat.py
+++ b/test/FIAT/unit/test_fiat.py
@@ -59,7 +59,7 @@
 from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1          # noqa: F401
 from FIAT.guzman_neilan import GuzmanNeilanSecondKindH1         # noqa: F401
 from FIAT.johnson_mercier import JohnsonMercier                 # noqa: F401
-from FIAT.bubble import Bubble
+from FIAT.bubble import Bubble, FacetBubble                     # noqa: F401
 from FIAT.enriched import EnrichedElement                       # noqa: F401
 from FIAT.nodal_enriched import NodalEnrichedElement
 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen    # noqa: F401
@@ -320,8 +320,15 @@ def __init__(self, a, b):
     "Histopolation(I, 1)",
     "Histopolation(I, 2)",
     "Bubble(I, 2)",
+    "Bubble(I, 2, 'integral')",
     "Bubble(T, 3)",
+    "Bubble(T, 3, 'integral')",
     "Bubble(S, 4)",
+    "Bubble(S, 4, 'integral')",
+    "FacetBubble(T, 2)",
+    "FacetBubble(T, 2, 'integral')",
+    "FacetBubble(S, 3)",
+    "FacetBubble(S, 3, 'integral')",
     "RestrictedElement(Lagrange(I, 2), restriction_domain='facet')",
     "RestrictedElement(Lagrange(T, 2), restriction_domain='vertex')",
     "RestrictedElement(Lagrange(T, 3), restriction_domain='facet')",
diff --git a/test/FIAT/unit/test_hdivtrace.py b/test/FIAT/unit/test_hdivtrace.py
index ab66f0369..ec405a190 100644
--- a/test/FIAT/unit/test_hdivtrace.py
+++ b/test/FIAT/unit/test_hdivtrace.py
@@ -25,7 +25,8 @@
 
 @pytest.mark.parametrize("dim", (2, 3))
 @pytest.mark.parametrize("degree", range(7))
-def test_basis_values(dim, degree):
+@pytest.mark.parametrize("variant", ("spectral", "integral"))
+def test_basis_values(dim, degree, variant):
     """Ensure that integrating simple monomials produces the expected results
     for each facet entity of the reference triangle and tetrahedron.
 
@@ -39,7 +40,7 @@ def test_basis_values(dim, degree):
 
     ref_el = ufc_simplex(dim)
     quadrule = make_quadrature(ufc_simplex(dim - 1), degree + 1)
-    fiat_element = HDivTrace(ref_el, degree)
+    fiat_element = HDivTrace(ref_el, degree, variant=variant)
     facet_element = fiat_element.dg_elements[dim - 1]
     nf = facet_element.space_dimension()
 

From e04a204843d7585154137aa0b0848e4066fcc908 Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Mon, 30 Jun 2025 14:17:44 +0100
Subject: [PATCH 05/11] cleanup

---
 FIAT/bubble.py     |  5 ++++-
 FIAT/hdiv_trace.py | 37 ++++++++++++++++++-------------------
 2 files changed, 22 insertions(+), 20 deletions(-)

diff --git a/FIAT/bubble.py b/FIAT/bubble.py
index 8d685b46e..563b3192d 100644
--- a/FIAT/bubble.py
+++ b/FIAT/bubble.py
@@ -16,7 +16,10 @@ class CodimBubble(RestrictedElement):
     """Bubbles of a certain codimension."""
 
     def __init__(self, ref_el, degree, codim, variant=None):
-        CG = IntegratedLegendre if variant == "integral" else Lagrange
+        if variant and variant.startswith("integral"):
+            CG = IntegratedLegendre
+        else:
+            CG = Lagrange
         element = CG(ref_el, degree, variant=variant)
 
         cell_dim = ref_el.get_dimension()
diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py
index 54b008844..ea7870243 100644
--- a/FIAT/hdiv_trace.py
+++ b/FIAT/hdiv_trace.py
@@ -81,29 +81,26 @@ def __init__(self, ref_el, degree, variant=None):
         topology = ref_el.get_topology()
         entity_dofs = {dim: {entity: [] for entity in topology[dim]} for dim in topology}
 
-        as_int = lambda x: sum(x) if isinstance(x, tuple) else x
-        # construct the DG element for the facets
+        # Construct the DG element for the facets
         dg_elements = {}
         for dim in topology:
-            if as_int(dim) == facet_sd:
+            fdim = sum(dim) if isinstance(dim, tuple) else dim
+            if fdim == facet_sd:
                 cell = ref_el.construct_subelement(dim)
                 dg_elements[dim] = construct_dg_element(cell, degree, variant)
 
         # Compute the dof numbering for all facet entities
         # and extract nodes
-        offset = 0
         nodes = []
         for facet_dim in sorted(dg_elements):
             element = dg_elements[facet_dim]
-            nf = element.space_dimension()
-            num_facets = len(topology[facet_dim])
-
-            for i in range(num_facets):
-                entity_dofs[facet_dim][i] = list(range(offset, offset + nf))
-                offset += nf
-                nodes.extend(transform_nodes(element.dual_basis(), ref_el, facet_dim, i))
+            facet_nodes = element.dual_basis()
+            for i in sorted(topology[facet_dim]):
+                cur = len(nodes)
+                nodes.extend(transform_nodes(facet_nodes, ref_el, facet_dim, i))
+                entity_dofs[facet_dim][i] = list(range(cur, len(nodes)))
 
-        # Setting up dual basis - only point evaluations
+        # Setting up dual basis
         dual = DualSet(nodes, ref_el, entity_dofs)
 
         # Degree of the element
@@ -193,13 +190,13 @@ def tabulate(self, order, points, entity=None):
                         phivals[key] = TraceError(msg)
 
             # Otherwise, extract non-zero values and insertion indices
+            element = self.dg_elements[facet_sd]
+            nf = element.space_dimension()
             for facet, ipts in facet_to_pts.items():
                 # Map points to the reference facet
                 new_points = map_to_reference_facet(points[ipts], vertices, facet)
 
                 # Retrieve values by tabulating the DG element
-                element = self.dg_elements[facet_sd]
-                nf = element.space_dimension()
                 nonzerovals = element.tabulate(order, new_points)[(0,)*facet_sd]
                 indices = slice(nf * facet, nf * (facet + 1))
 
@@ -220,12 +217,11 @@ def tabulate(self, order, points, entity=None):
                 # Retrieve function evaluations (order = 0 case)
                 offset = 0
                 for facet_dim in sorted(self.dg_elements):
-                    element = self.dg_elements[facet_dim]
-                    nf = element.space_dimension()
-
                     # Loop over the number of facets until we find a facet
                     # with matching dimension and id
-                    for i in self.ref_el.get_topology()[facet_dim]:
+                    element = self.dg_elements[facet_dim]
+                    nf = element.space_dimension()
+                    for i in sorted(self.ref_el.get_topology()[facet_dim]):
                         # Found it! Grab insertion indices
                         if (facet_dim, i) == entity:
                             nonzerovals = element.tabulate(0, points)[(0,)*facet_sd]
@@ -259,7 +255,10 @@ def construct_dg_element(ref_el, degree, variant):
     """Constructs a discontinuous galerkin element of a given degree
     on a particular reference cell.
     """
-    DG = Legendre if variant == "integral" else DiscontinuousLagrange
+    if variant and variant.startswith("integral"):
+        DG = Legendre
+    else:
+        DG = DiscontinuousLagrange
     if ref_el.get_shape() in [LINE, TRIANGLE]:
         dg_element = DG(ref_el, degree, variant)
 

From a129b96be3093c7d87d1716cba99dfe13ec79102 Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Mon, 30 Jun 2025 20:57:49 +0100
Subject: [PATCH 06/11] Implement HDivTrace as a macroelement

---
 FIAT/bubble.py                   |   7 +-
 FIAT/discontinuous_lagrange.py   |   5 +-
 FIAT/dual_set.py                 |   8 +-
 FIAT/expansions.py               |  62 ++++-
 FIAT/hdiv_trace.py               | 436 ++++++-------------------------
 FIAT/hierarchical.py             |   7 +-
 FIAT/lagrange.py                 |   6 +
 FIAT/polynomial_set.py           |   7 +-
 FIAT/reference_element.py        |  61 ++++-
 finat/element_factory.py         |  27 +-
 finat/ufl/finiteelement.py       |  15 ++
 test/FIAT/unit/test_fiat.py      |  34 ++-
 test/FIAT/unit/test_hdivtrace.py |  31 ++-
 13 files changed, 264 insertions(+), 442 deletions(-)

diff --git a/FIAT/bubble.py b/FIAT/bubble.py
index 563b3192d..6ecfeb826 100644
--- a/FIAT/bubble.py
+++ b/FIAT/bubble.py
@@ -7,7 +7,6 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 from FIAT.lagrange import Lagrange
-from FIAT.hierarchical import IntegratedLegendre
 from FIAT.restricted import RestrictedElement
 from itertools import chain
 
@@ -16,11 +15,7 @@ class CodimBubble(RestrictedElement):
     """Bubbles of a certain codimension."""
 
     def __init__(self, ref_el, degree, codim, variant=None):
-        if variant and variant.startswith("integral"):
-            CG = IntegratedLegendre
-        else:
-            CG = Lagrange
-        element = CG(ref_el, degree, variant=variant)
+        element = Lagrange(ref_el, degree, variant=variant)
 
         cell_dim = ref_el.get_dimension()
         assert cell_dim == max(element.entity_dofs().keys())
diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py
index fd08f8fa3..db5c06b83 100644
--- a/FIAT/discontinuous_lagrange.py
+++ b/FIAT/discontinuous_lagrange.py
@@ -14,6 +14,7 @@
 from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points
 from FIAT.polynomial_set import mis
 from FIAT.check_format_variant import parse_lagrange_variant
+from FIAT.hierarchical import Legendre
 
 
 def make_entity_permutations(dim, npoints):
@@ -215,6 +216,8 @@ class DiscontinuousLagrange(finite_element.CiarletElement):
                               macroelement for Scott-Vogelius.
     """
     def __new__(cls, ref_el, degree, variant="equispaced"):
+        if variant and variant.startswith("integral"):
+            return Legendre(ref_el, degree, variant=variant)
         if degree == 0:
             splitting, _ = parse_lagrange_variant(variant, discontinuous=True)
             if splitting is None and not ref_el.is_macrocell():
@@ -230,7 +233,7 @@ def __init__(self, ref_el, degree, variant="equispaced"):
             dual = BrokenLagrangeDualSet(ref_el, degree, point_variant=point_variant)
         else:
             dual = DiscontinuousLagrangeDualSet(ref_el, degree, point_variant=point_variant)
-        if ref_el.shape == LINE:
+        if ref_el.shape == LINE and not ref_el.is_trace():
             # In 1D we can use the primal basis as the expansion set,
             # avoiding any round-off coming from a basis transformation
             points = get_lagrange_points(dual)
diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py
index c09db974f..539c87553 100644
--- a/FIAT/dual_set.py
+++ b/FIAT/dual_set.py
@@ -277,15 +277,15 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations):
     parent_cell = ref_el.get_parent()
     if parent_cell is None:
         return nodes, ref_el, entity_ids, entity_permutations
-    parent_ids = {}
+    parent_top = parent_cell.get_topology()
+    parent_ids = {dim: {entity: [] for entity in parent_top[dim]} for dim in parent_top}
     parent_permutations = None
     parent_to_children = ref_el.get_parent_to_children()
 
-    if all(isinstance(node, functional.PointEvaluation) for node in nodes):
+    if all(isinstance(node, functional.PointEvaluation) for node in nodes) and not ref_el.is_trace():
         # Merge Lagrange dual with lexicographical reordering
         parent_nodes = []
         for dim in sorted(parent_to_children):
-            parent_ids[dim] = {}
             for entity in sorted(parent_to_children[dim]):
                 cur = len(parent_nodes)
                 for child_dim, child_entity in parent_to_children[dim][entity]:
@@ -296,9 +296,7 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations):
         # Merge everything else with the same node ordering
         parent_nodes = nodes
         for dim in sorted(parent_to_children):
-            parent_ids[dim] = {}
             for entity in sorted(parent_to_children[dim]):
-                parent_ids[dim][entity] = []
                 for child_dim, child_entity in parent_to_children[dim][entity]:
                     parent_ids[dim][entity].extend(entity_ids[child_dim][child_entity])
 
diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index 69d0d4b5a..959490552 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -268,7 +268,7 @@ def __init__(self, ref_el, scale=None, variant=None):
         base_ref_el = reference_element.default_simplex(sd)
         base_verts = base_ref_el.get_vertices()
 
-        self.affine_mappings = [reference_element.make_affine_mapping(
+        self.affine_mappings = [get_affine_mapping(
                                 ref_el.get_vertices_of_subcomplex(top[sd][cell]),
                                 base_verts) for cell in top[sd]]
         if scale is None:
@@ -345,6 +345,8 @@ def distance(alpha, beta):
 
     def _tabulate(self, n, pts, order=0):
         """A version of tabulate() that also works for a single point."""
+        from FIAT.hdiv_trace import TraceError
+        from FIAT.polynomial_set import mis
         pts = numpy.asarray(pts)
         unique = self.continuity is not None and order == 0
         cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=unique)
@@ -354,6 +356,28 @@ def _tabulate(self, n, pts, order=0):
         if not self.ref_el.is_macrocell():
             return phis[0]
 
+        if self.ref_el.is_trace():
+            parent = self.ref_el.get_parent()
+            tdim = self.ref_el.get_spatial_dimension()
+            gdim = parent.get_spatial_dimension()
+            for cell in phis:
+                # Promote facet keys to cell keys
+                phi = phis[cell][(0,) * tdim]
+                # Raise TraceError on gradient tabulations
+                msg = "Gradients on trace elements are not well-defined."
+                phis[cell] = {alpha: phi if sum(alpha) == 0 else TraceError(msg)
+                              for i in range(order+1)
+                              for alpha in mis(gdim, i)}
+            if tdim == 0 and len(phis) == 0:
+                # Hack for TensorProduct HDivTrace: do not raise TraceError on the interval
+                for cell in parent.topology[gdim]:
+                    phis[cell] = {(0,)*gdim: numpy.zeros((1, 1))}
+            elif sum(len(cell_point_map[cell]) for cell in cell_point_map) < len(pts):
+                # Raise TraceError when interior points fail to be binned on facets
+                for cell in parent.topology[gdim]:
+                    msg = "The HDivTrace element can only be tabulated on facets."
+                    phis[cell] = {(0,)*gdim: TraceError(msg)}
+
         if pts.dtype == object:
             # If binning is undefined, scale by the characteristic function of each subcell
             Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique)
@@ -362,13 +386,14 @@ def _tabulate(self, n, pts, order=0):
                     phi[alpha] *= Xi[cell]
         elif not unique:
             # If binning is not unique, divide by the multiplicity of each point
-            mult = numpy.zeros(pts.shape[:-1])
+            mult = numpy.zeros(pts.shape[:-1], dtype=int)
             for cell, ipts in cell_point_map.items():
                 mult[ipts] += 1
-            for cell, ipts in cell_point_map.items():
-                phi = phis[cell]
-                for alpha in phi:
-                    phi[alpha] /= mult[None, ipts]
+            if (mult != 1).any():
+                for cell, ipts in cell_point_map.items():
+                    phi = phis[cell]
+                    for alpha in phi:
+                        phi[alpha] /= mult[None, ipts]
 
         # Insert subcell tabulations into the corresponding submatrices
         idx = lambda *args: args if args[1] is Ellipsis else numpy.ix_(*args)
@@ -377,6 +402,9 @@ def _tabulate(self, n, pts, order=0):
         result = {}
         base_phi = tuple(phis.values())[0]
         for alpha in base_phi:
+            if isinstance(base_phi[alpha], TraceError):
+                result[alpha] = base_phi[alpha]
+                continue
             dtype = base_phi[alpha].dtype
             result[alpha] = numpy.zeros((num_phis, *pts.shape[:-1]), dtype=dtype)
             for cell in cell_point_map:
@@ -497,8 +525,7 @@ def get_dmats(self, degree, cell=0):
     def tabulate(self, n, pts):
         if len(pts) == 0:
             return numpy.array([])
-        sd = self.ref_el.get_spatial_dimension()
-        return self._tabulate(n, pts)[(0,) * sd]
+        return tuple(self._tabulate(n, pts).values())[0]
 
     def tabulate_derivatives(self, n, pts):
         from FIAT.polynomial_set import mis
@@ -591,13 +618,30 @@ def __init__(self, ref_el, **kwargs):
         super().__init__(ref_el, **kwargs)
 
 
+def get_affine_mapping(xs, ys):
+    """Constructs (A,b) such that x --> A * x + b is the affine
+    mapping from the simplex defined by xs to the simplex defined by ys.
+    Uses least-squares when the simplex dimension does not match the spatial
+    dimension.
+    """
+    X = numpy.asarray(xs)
+    Y = numpy.asarray(ys)
+    A = X[1:] - X[:1]
+    B = Y[1:] - Y[:1]
+    if A.shape[0] == A.shape[1]:
+        C = numpy.linalg.solve(A, B)
+    else:
+        C, *_ = numpy.linalg.lstsq(A, B)
+    b = Y[0] - numpy.dot(X[0], C)
+    return C.T, b
+
+
 def polynomial_dimension(ref_el, n, continuity=None):
     """Returns the dimension of the space of polynomials of degree no
     greater than n on the reference complex."""
     if ref_el.get_shape() == reference_element.POINT:
         if n > 0:
             raise ValueError("Only degree zero polynomials supported on point elements.")
-        return 1
     top = ref_el.get_topology()
     if continuity == "C0":
         space_dimension = sum(math.comb(n - 1, dim) * len(top[dim]) for dim in top)
diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py
index ea7870243..2ea70a4ee 100644
--- a/FIAT/hdiv_trace.py
+++ b/FIAT/hdiv_trace.py
@@ -4,24 +4,79 @@
 #
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
-import numpy as np
-from collections import defaultdict
-from FIAT.barycentric_interpolation import get_lagrange_points
 from FIAT.discontinuous_lagrange import DiscontinuousLagrange
-from FIAT.hierarchical import Legendre
-from FIAT.dual_set import DualSet
-from FIAT.finite_element import FiniteElement
-from FIAT.functional import IntegralMoment, PointEvaluation
-from FIAT.polynomial_set import mis
-from FIAT.quadrature import FacetQuadratureRule
-from FIAT.reference_element import (ufc_simplex, POINT,
-                                    LINE, QUADRILATERAL,
-                                    TRIANGLE, TETRAHEDRON,
-                                    TENSORPRODUCT)
-from FIAT.tensor_product import TensorProductElement
+from FIAT.reference_element import SimplicialComplex
+
+
+class TraceSimplicialComplex(SimplicialComplex):
+
+    def __init__(self, parent, codim=1):
+        sd = parent.get_spatial_dimension() - codim
+        parent_top = parent.topology
+        topology = dict(parent_top)
+        for dim in parent_top:
+            if dim > sd:
+                topology.pop(dim)
+
+        facet = parent.construct_subelement(sd)
+        self.base_ref_el = facet
+        self._parent_complex = parent
+        self._parent_simplex = parent
+
+        # dict mapping child facets to their parent facet
+        parent_to_children = {d: {e: [(d, e)]
+                              for e in topology[d]}
+                              for d in topology}
+        # dict mapping parent facets to their children
+        child_to_parent = {d: {e: (d, e)
+                           for e in topology[d]}
+                           for d in topology}
+        self._child_to_parent = child_to_parent
+        self._parent_to_children = parent_to_children
+
+        # dict mapping cells to their boundary facets for each dimension,
+        # while respecting the ordering on the parent simplex
+        connectivity = {cell: {dim: [] for dim in topology} for cell in topology[sd]}
+        for cell in topology[sd]:
+            for dim in topology:
+                for entity in topology[dim]:
+                    if set(topology[dim][entity]) <= set(topology[sd][cell]):
+                        connectivity[cell][dim].append(entity)
+        self._cell_connectivity = connectivity
+
+        # dict mapping subentity dimension to interior facets
+        interior_facets = {dim: [entity for entity in child_to_parent[dim]
+                                 if child_to_parent[dim][entity][0] == sd]
+                           for dim in sorted(child_to_parent)}
+        self._interior_facets = interior_facets
+        super().__init__(facet.shape, parent.vertices, topology)
+
+    def get_spatial_dimension(self):
+        return self.base_ref_el.get_spatial_dimension()
+
+    def get_dimension(self):
+        return self.base_ref_el.get_dimension()
+
+    def construct_subelement(self, dim):
+        return self._parent_simplex.construct_subelement(dim)
+
+    def get_parent(self):
+        return self._parent_simplex
+
+    def get_parent_complex(self):
+        return self._parent_complex
+
+    def get_parent_to_children(self):
+        return self._parent_to_children
+
+    def get_cell_connectivity(self):
+        return self._cell_connectivity
+
+    def is_trace(self):
+        return True
 
-# Numerical tolerance for facet-entity identifications
-epsilon = 1e-10
+    def is_macrocell(self):
+        return True
 
 
 class TraceError(Exception):
@@ -33,7 +88,7 @@ def __init__(self, msg):
         self.msg = msg
 
 
-class HDivTrace(FiniteElement):
+class HDivTrace(DiscontinuousLagrange):
     """Class implementing the trace of hdiv elements. This class
     is a stand-alone element family that produces a DG-facet field.
     This element is what's produced after performing the trace
@@ -42,8 +97,7 @@ class HDivTrace(FiniteElement):
     This element is also known as the discontinuous trace field that
     arises in several DG formulations.
     """
-
-    def __init__(self, ref_el, degree, variant=None):
+    def __new__(cls, ref_el, degree, variant="equispaced_interior"):
         """Constructor for the HDivTrace element.
 
         :arg ref_el: A reference element, which may be a tensor product
@@ -53,345 +107,5 @@ def __init__(self, ref_el, degree, variant=None):
                      varying degrees.
         :arg variant: The point distribution variant passed on to recursivenodes.
         """
-        sd = ref_el.get_spatial_dimension()
-        if sd in (0, 1):
-            raise ValueError("Cannot take the trace of a %d-dim cell." % sd)
-
-        # Store the degrees if on a tensor product cell
-        if ref_el.get_shape() == TENSORPRODUCT:
-            try:
-                degree = tuple(degree)
-            except TypeError:
-                degree = (degree,) * len(ref_el.cells)
-
-            assert len(ref_el.cells) == len(degree), (
-                "Number of specified degrees must be equal to the number of cells."
-            )
-        else:
-            if ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON, QUADRILATERAL]:
-                raise NotImplementedError(
-                    "Trace element on a %s not implemented" % type(ref_el)
-                )
-            # Cannot have varying degrees for these reference cells
-            if isinstance(degree, tuple):
-                raise ValueError("Must have a tensor product cell if providing multiple degrees")
-
-        # Initialize entity dofs
-        facet_sd = sd - 1
-        topology = ref_el.get_topology()
-        entity_dofs = {dim: {entity: [] for entity in topology[dim]} for dim in topology}
-
-        # Construct the DG element for the facets
-        dg_elements = {}
-        for dim in topology:
-            fdim = sum(dim) if isinstance(dim, tuple) else dim
-            if fdim == facet_sd:
-                cell = ref_el.construct_subelement(dim)
-                dg_elements[dim] = construct_dg_element(cell, degree, variant)
-
-        # Compute the dof numbering for all facet entities
-        # and extract nodes
-        nodes = []
-        for facet_dim in sorted(dg_elements):
-            element = dg_elements[facet_dim]
-            facet_nodes = element.dual_basis()
-            for i in sorted(topology[facet_dim]):
-                cur = len(nodes)
-                nodes.extend(transform_nodes(facet_nodes, ref_el, facet_dim, i))
-                entity_dofs[facet_dim][i] = list(range(cur, len(nodes)))
-
-        # Setting up dual basis
-        dual = DualSet(nodes, ref_el, entity_dofs)
-
-        # Degree of the element
-        deg = max(e.degree() for e in dg_elements.values())
-
-        super().__init__(ref_el, dual, order=deg,
-                         formdegree=facet_sd,
-                         mapping="affine")
-
-        # Set up facet elements
-        self.dg_elements = dg_elements
-
-        # Degree for quadrature rule
-        self.polydegree = deg
-
-    def degree(self):
-        """Return the degree of the (embedding) polynomial space."""
-        return self.polydegree
-
-    def get_nodal_basis(self):
-        """Return the nodal basis, encoded as a PolynomialSet object,
-        for the finite element."""
-        raise NotImplementedError("get_nodal_basis not implemented for the trace element.")
-
-    def get_coeffs(self):
-        """Return the expansion coefficients for the basis of the
-        finite element."""
-        raise NotImplementedError("get_coeffs not implemented for the trace element.")
-
-    def tabulate(self, order, points, entity=None):
-        """Return tabulated values of derivatives up to a given order of
-        basis functions at given points.
-
-        :arg order: The maximum order of derivative.
-        :arg points: An iterable of points.
-        :arg entity: Optional (dimension, entity number) pair
-                     indicating which topological entity of the
-                     reference element to tabulate on.  If ``None``,
-                     tabulated values are computed by geometrically
-                     approximating which facet the points are on.
-
-        .. note ::
-
-           Performing illegal tabulations on this element will result in either
-           a tabulation table of `numpy.nan` arrays (`entity=None` case), or
-           insertions of the `TraceError` exception class. This is due to the
-           fact that performing cell-wise tabulations, or asking for any order
-           of derivative evaluations, are not mathematically well-defined.
-        """
-        sd = self.ref_el.get_spatial_dimension()
-        facet_sd = sd - 1
-        evalkey = (0,) * sd
-
-        # Initializing dictionary with zeros
-        phivals = {}
-        for i in range(order + 1):
-            alphas = mis(sd, i)
-            for alpha in alphas:
-                phivals[alpha] = np.zeros(shape=(self.space_dimension(), len(points)))
-                if alpha != evalkey:
-                    # If asking for gradient evaluations, insert TraceError in gradient slots
-                    phivals[alpha] = TraceError("Gradients on trace elements are not well-defined.")
-
-        # If entity is None, identify facet using numerical tolerance and
-        # return the tabulated values
-        if entity is None or entity == (sd, 0):
-            # NOTE: Numerical approximation of the facet id is currently only
-            # implemented for simplex reference cells.
-            if self.ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON]:
-                raise NotImplementedError(
-                    "Tabulating this element on a %s cell without providing "
-                    "an entity is not currently supported." % type(self.ref_el)
-                )
-
-            # Attempt to identify which facet (if any) the given points are on
-            vertices = self.ref_el.vertices
-            coordinates = barycentric_coordinates(points, vertices)
-            facet_to_pts, success = extract_facets(coordinates)
-
-            # If not successful, return NaNs
-            if not success:
-                for key in phivals:
-                    if entity is None:
-                        phivals[key].fill(np.nan)
-                    else:
-                        msg = "The HDivTrace element can only be tabulated on facets."
-                        phivals[key] = TraceError(msg)
-
-            # Otherwise, extract non-zero values and insertion indices
-            element = self.dg_elements[facet_sd]
-            nf = element.space_dimension()
-            for facet, ipts in facet_to_pts.items():
-                # Map points to the reference facet
-                new_points = map_to_reference_facet(points[ipts], vertices, facet)
-
-                # Retrieve values by tabulating the DG element
-                nonzerovals = element.tabulate(order, new_points)[(0,)*facet_sd]
-                indices = slice(nf * facet, nf * (facet + 1))
-
-                # Insert non-zero values in appropriate place
-                phivals[evalkey][indices, ipts] = nonzerovals
-
-        else:
-            entity_dim, _ = entity
-
-            # If the user is directly specifying cell-wise tabulation, return
-            # TraceErrors in dict for appropriate handling in the form compiler
-            if entity_dim not in self.dg_elements:
-                for key in phivals:
-                    msg = "The HDivTrace element can only be tabulated on facets."
-                    phivals[key] = TraceError(msg)
-
-            else:
-                # Retrieve function evaluations (order = 0 case)
-                offset = 0
-                for facet_dim in sorted(self.dg_elements):
-                    # Loop over the number of facets until we find a facet
-                    # with matching dimension and id
-                    element = self.dg_elements[facet_dim]
-                    nf = element.space_dimension()
-                    for i in sorted(self.ref_el.get_topology()[facet_dim]):
-                        # Found it! Grab insertion indices
-                        if (facet_dim, i) == entity:
-                            nonzerovals = element.tabulate(0, points)[(0,)*facet_sd]
-                            indices = slice(offset, offset + nf)
-                        offset += nf
-
-                # Insert non-zero values in appropriate place
-                phivals[evalkey][indices] = nonzerovals
-
-        return phivals
-
-    def value_shape(self):
-        """Return the value shape of the finite element functions."""
-        return ()
-
-    def dmats(self):
-        """Return dmats: expansion coefficients for basis function
-        derivatives."""
-        raise NotImplementedError("dmats not implemented for the trace element.")
-
-    def get_num_members(self, arg):
-        """Return number of members of the expansion set."""
-        raise NotImplementedError("get_num_members not implemented for the trace element.")
-
-    @staticmethod
-    def is_nodal():
-        return True
-
-
-def construct_dg_element(ref_el, degree, variant):
-    """Constructs a discontinuous galerkin element of a given degree
-    on a particular reference cell.
-    """
-    if variant and variant.startswith("integral"):
-        DG = Legendre
-    else:
-        DG = DiscontinuousLagrange
-    if ref_el.get_shape() in [LINE, TRIANGLE]:
-        dg_element = DG(ref_el, degree, variant)
-
-    # Quadrilateral facets could be on a FiredrakeQuadrilateral.
-    # In this case, we treat this as an interval x interval cell:
-    elif ref_el.get_shape() == QUADRILATERAL:
-        dg_line = DG(ufc_simplex(1), degree, variant)
-        dg_element = TensorProductElement(dg_line, dg_line)
-
-    # This handles the more general case for facets:
-    elif ref_el.get_shape() == TENSORPRODUCT:
-        assert len(degree) == len(ref_el.cells), (
-            "Must provide the same number of degrees as the number "
-            "of cells that make up the tensor product cell."
-        )
-        sub_elements = [construct_dg_element(c, d, variant)
-                        for c, d in zip(ref_el.cells, degree)
-                        if c.get_shape() != POINT]
-
-        if len(sub_elements) > 1:
-            dg_element = TensorProductElement(*sub_elements)
-        else:
-            dg_element, = sub_elements
-
-    else:
-        raise NotImplementedError(
-            "Reference cells of type %s not currently supported" % type(ref_el)
-        )
-
-    return dg_element
-
-
-def transform_nodes(ells, ref_el, facet_dim, facet_id):
-    """Map functionals into a given facet."""
-    try:
-        facet_pts = get_lagrange_points(ells)
-        transform = ref_el.get_entity_transform(facet_dim, facet_id)
-        pts = transform(facet_pts)
-        for pt in pts:
-            yield PointEvaluation(ref_el, pt)
-    except ValueError:
-        Q_ref, = set(ell.Q for ell in ells)
-        Q = FacetQuadratureRule(ref_el, facet_dim, facet_id, Q_ref)
-        for ell in ells:
-            yield IntegralMoment(ref_el, Q, ell.f_at_qpts)
-
-
-# The following functions are credited to Marie E. Rognes:
-def extract_facets(coordinates, tolerance=epsilon):
-    """Determines whether a set of points (described in barycentric coordinates)
-    are all on facet sub-entities, and return a dict mapping facets to
-    point indices and whether the search has been successful.
-
-    :arg coordinates: A set of points described in barycentric coordinates.
-    :arg tolerance: A fixed tolerance for geometric identifications.
-    """
-    facet_to_pts = defaultdict(list)
-    for ipt, c in enumerate(coordinates):
-        on_facet = set(i for (i, l) in enumerate(c) if abs(l) < tolerance)
-        try:
-            f, = on_facet
-        except ValueError:
-            # Handle coordinates not on facets
-            return ({}, False)
-        facet_to_pts[f].append(ipt)
-
-    # If all points are on facets, return indices and success
-    return (facet_to_pts, True)
-
-
-def barycentric_coordinates(points, vertices):
-    """Computes the barycentric coordinates for a set of points relative to a
-    simplex defined by a set of vertices.
-
-    :arg points: A set of points.
-    :arg vertices: A set of vertices that define the simplex.
-    """
-
-    # Form mapping matrix
-    T = (np.asarray(vertices[:-1]) - vertices[-1]).T
-    invT = np.linalg.inv(T)
-
-    points = np.asarray(points)
-    bary = np.einsum("ij,kj->ki", invT, (points - vertices[-1]))
-    last = (1 - bary.sum(axis=1))
-    return np.concatenate([bary, last[..., np.newaxis]], axis=1)
-
-
-def map_from_reference_facet(point, vertices):
-    """Evaluates the physical coordinate of a point using barycentric
-    coordinates.
-
-    :arg point: The reference points to be mapped to the facet.
-    :arg vertices: The vertices defining the physical element.
-    """
-
-    # Compute the barycentric coordinates of the point relative to the reference facet
-    reference_simplex = ufc_simplex(len(vertices) - 1)
-    reference_vertices = reference_simplex.get_vertices()
-    coords = barycentric_coordinates([point, ], reference_vertices)[0]
-
-    # Evaluates the physical coordinate of the point using barycentric coordinates
-    point = sum(vertices[j] * coords[j] for j in range(len(coords)))
-    return tuple(point)
-
-
-def map_to_reference_facet(points, vertices, facet):
-    """Given a set of points and vertices describing a facet of a simplex in n-dimensional
-    coordinates (where the points lie on the facet), map the points to the reference simplex
-    of dimension (n-1).
-
-    :arg points: A set of points in n-D.
-    :arg vertices: A set of vertices describing a facet of a simplex in n-D.
-    :arg facet: Integer representing the facet number.
-    """
-
-    # Compute the barycentric coordinates of the points with respect to the
-    # full physical simplex
-    all_coords = barycentric_coordinates(points, vertices)
-
-    # Extract vertices of the reference facet
-    reference_facet_simplex = ufc_simplex(len(vertices) - 2)
-    reference_vertices = reference_facet_simplex.get_vertices()
-
-    reference_points = []
-    for (i, coords) in enumerate(all_coords):
-        # Extract the correct subset of barycentric coordinates since we know
-        # which facet we are on
-        new_coords = [coords[j] for j in range(len(coords)) if j != facet]
-
-        # Evaluate the reference coordinate of a point in barycentric coordinates
-        reference_pt = sum(np.asarray(reference_vertices[j]) * new_coords[j]
-                           for j in range(len(new_coords)))
-
-        reference_points += [reference_pt]
-    return reference_points
+        facets = TraceSimplicialComplex(ref_el)
+        return DiscontinuousLagrange(facets, degree, variant=variant)
diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py
index 17c415710..535ea7d7b 100644
--- a/FIAT/hierarchical.py
+++ b/FIAT/hierarchical.py
@@ -51,13 +51,14 @@ def __init__(self, ref_el, degree, codim=0):
                 continue
 
             ref_facet = ref_el.construct_subelement(dim)
-            poly_set = ONPolynomialSet(ref_facet, degree)
+            poly_set = ONPolynomialSet(ref_facet, degree, scale=1)
             Q_ref = create_quadrature(ref_facet, 2 * degree)
             phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim]
             for entity in sorted(top[dim]):
                 cur = len(nodes)
+                Jdet = ref_el.volume_of_subcomplex(dim, entity)
                 Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-                nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
+                nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi/Jdet) for phi in phis)
                 entity_ids[dim][entity] = list(range(cur, len(nodes)))
                 entity_permutations[dim][entity] = perms
 
@@ -69,7 +70,7 @@ class Legendre(finite_element.CiarletElement):
     def __new__(cls, ref_el, degree, variant=None):
         if degree == 0:
             splitting, _ = parse_lagrange_variant(variant, integral=True)
-            if splitting is None:
+            if splitting is None and not ref_el.is_macrocell():
                 # FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size()
                 return P0.P0(ref_el)
         return super().__new__(cls)
diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py
index eead0cfa6..d8219a4ce 100644
--- a/FIAT/lagrange.py
+++ b/FIAT/lagrange.py
@@ -10,6 +10,7 @@
 from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points
 from FIAT.reference_element import LINE
 from FIAT.check_format_variant import parse_lagrange_variant
+from FIAT.hierarchical import IntegratedLegendre
 
 
 class LagrangeDualSet(dual_set.DualSet):
@@ -72,6 +73,11 @@ class Lagrange(finite_element.CiarletElement):
                         entity id. The DOFs are always sorted by the entity ordering
                         and then lexicographically by lattice multiindex.
     """
+    def __new__(cls, ref_el, degree, variant="equispaced", sort_entities=False):
+        if variant and variant.startswith("integral"):
+            return IntegratedLegendre(ref_el, degree, variant=variant)
+        return super().__new__(cls)
+
     def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False):
         splitting, point_variant = parse_lagrange_variant(variant)
         if splitting is not None:
diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py
index 534c18642..2d158fce4 100644
--- a/FIAT/polynomial_set.py
+++ b/FIAT/polynomial_set.py
@@ -68,8 +68,11 @@ def tabulate_new(self, pts):
 
     def tabulate(self, pts, jet_order=0):
         """Returns the values of the polynomial set."""
-        base_vals = self.expansion_set._tabulate(self.embedded_degree, pts, order=jet_order)
-        result = {alpha: numpy.dot(self.coeffs, base_vals[alpha]) for alpha in base_vals}
+        result = self.expansion_set._tabulate(self.embedded_degree, pts, order=jet_order)
+        for alpha in result:
+            if isinstance(result[alpha], Exception):
+                continue
+            result[alpha] = numpy.dot(self.coeffs, result[alpha])
         return result
 
     def get_expansion_set(self):
diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index 547da2c04..8a2aed15e 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -94,6 +94,8 @@ def make_lattice(verts, n, interior=0, variant=None):
     family = _decode_family(family)
     D = len(verts)
     X = numpy.array(verts)
+    if D == 1 and n == 0:
+        return list(map(tuple, X))
     get_point = lambda alpha: tuple(numpy.dot(_recursive(D - 1, n, alpha, family), X))
     return list(map(get_point, multiindex_equal(D, n, interior)))
 
@@ -323,6 +325,9 @@ def extrinsic_orientation_permutation_map(self):
     def is_simplex(self):
         return False
 
+    def is_trace(self):
+        return False
+
     def is_macrocell(self):
         return False
 
@@ -341,6 +346,10 @@ def get_parent_complex(self):
     def is_parent(self, other, strict=False):
         """Return whether this cell is the parent of the other cell."""
         parent = other
+        while self.is_trace():
+            self = self.get_parent_complex()
+        while parent.is_trace():
+            parent = parent.get_parent_complex()
         if strict:
             parent = parent.get_parent_complex()
         while parent is not None:
@@ -357,6 +366,8 @@ def __eq__(self, other):
             return False
         atop = self.get_topology()
         btop = other.get_topology()
+        if max(atop) != max(btop):
+            return False
         for dim in atop:
             if set(atop[dim].values()) != set(btop[dim].values()):
                 return False
@@ -566,7 +577,6 @@ def get_entity_transform(self, dim, entity):
         """
         topology = self.get_topology()
         celldim = self.get_spatial_dimension()
-        codim = celldim - dim
         if dim == 0:
             # Special case vertices.
             i, = topology[dim][entity]
@@ -578,7 +588,7 @@ def get_entity_transform(self, dim, entity):
         else:
             subcell = self.construct_subelement(dim)
             subdim = subcell.get_spatial_dimension()
-            assert subdim == celldim - codim
+            assert subdim == dim
 
             # Entity vertices in entity space.
             v_e = numpy.asarray(subcell.get_vertices())
@@ -610,19 +620,31 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
         if entity is None:
             entity = (self.get_spatial_dimension(), 0)
         entity_dim, entity_id = entity
-        top = self.get_topology()
-        sd = self.get_spatial_dimension()
+        sd = len(self.vertices[0])
 
         # get a subcell containing the entity and the restriction indices of the entity
+        outside = None
         indices = slice(None)
+        top = self.get_topology()
         subcomplex = top[entity_dim][entity_id]
         if entity_dim != sd:
-            cell_id = self.connectivity[(entity_dim, sd)][0][0]
-            indices = [i for i, v in enumerate(top[sd][cell_id]) if v in subcomplex]
-            subcomplex = top[sd][cell_id]
+            parent = self.get_parent() if self.is_trace() else self
+            cell_id = parent.connectivity[(entity_dim, sd)][0][0]
+            top = parent.get_topology()
+            subcell = top[sd][cell_id]
+            while len(subcell) > sd + 1:
+                # construct a simplex if we have a hypercube
+                k = max(set(subcell) - set(subcomplex))
+                subcell = subcell[:k] + subcell[k+1:]
+
+            indices = [i for i, v in enumerate(subcell) if v in subcomplex]
+            subcomplex = subcell
+
+            outside = [i for i in range(len(subcomplex)) if i not in indices]
+            indices.extend(outside)
 
         cell_verts = self.get_vertices_of_subcomplex(subcomplex)
-        ref_verts = numpy.eye(sd + 1)
+        ref_verts = numpy.eye(len(cell_verts))
         A, b = make_affine_mapping(cell_verts, ref_verts)
         A, b = A[indices], b[indices]
         if rescale:
@@ -630,8 +652,12 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
             h = 1 / numpy.linalg.norm(A, axis=1)
             b *= h
             A *= h[:, None]
-        out = numpy.dot(points, A.T)
-        return numpy.add(out, b, out=out)
+        bary = numpy.dot(points, A.T)
+        numpy.add(bary, b, out=bary)
+
+        if outside is not None:
+            bary[:, -len(outside):] = -abs(bary[:, -len(outside):])
+        return bary
 
     def compute_bubble(self, points, entity=None):
         """Returns the lowest-order bubble on an entity evaluated at the given
@@ -1249,9 +1275,9 @@ def distance_to_point_l1(self, point, rescale=False):
 
         For more information see the docstring for the UFCSimplex method."""
         subcell_dimensions = self.get_dimension()
-        assert len(point) == sum(subcell_dimensions)
-        point_slices = TensorProductCell._split_slices(subcell_dimensions)
         point = numpy.asarray(point)
+        assert point.shape[-1] == sum(subcell_dimensions)
+        point_slices = TensorProductCell._split_slices(subcell_dimensions)
         return sum(c.distance_to_point_l1(point[..., s], rescale=rescale)
                    for c, s in zip(self.cells, point_slices))
 
@@ -1390,6 +1416,17 @@ def construct_subelement(self, dimension):
             sub_element = self.product.construct_subelement((dimension,) + (0,)*(len(self.product.cells) - 1))
             return flatten_reference_cube(sub_element)
 
+    def get_facet_element(self):
+        sd = self.get_spatial_dimension()
+        return self.construct_subelement(sd-1)
+
+    def make_points(self, dim, entity, order, variant=None, interior=1):
+        if dim == 1:
+            entity_verts = self.get_vertices_of_subcomplex(self.topology[dim][entity])
+            return make_lattice(entity_verts, order, interior=interior, variant=variant)
+        else:
+            raise NotImplementedError("Cannot make points on non-simplex facets yet.")
+
     def get_entity_transform(self, dim, entity_i):
         """Returns a mapping of point coordinates from the
         `entity_i`-th subentity of dimension `dim` to the cell.
diff --git a/finat/element_factory.py b/finat/element_factory.py
index 51c4d3d66..fcca67d3e 100644
--- a/finat/element_factory.py
+++ b/finat/element_factory.py
@@ -158,19 +158,21 @@ def convert_finiteelement(element, **kwargs):
         codim = 1 if element.family() == "Boundary Quadrature" else 0
         return finat.make_quadrature_element(cell, degree, scheme, codim), set()
     lmbda = supported_elements[element.family()]
-    if element.family() == "Real" and element.cell.cellname() in {"quadrilateral", "hexahedron"}:
-        lmbda = None
-        element = finat.ufl.FiniteElement("DQ", element.cell, 0)
+    if element.cell.cellname() in supported_tensor_product_cells:
+        if element.family() == "Real":
+            lmbda = None
+            element = finat.ufl.FiniteElement("DQ", element.cell, 0)
+        elif element.family() == "HDiv Trace":
+            lmbda = None
+
     if lmbda is None:
-        if element.cell.cellname() == "quadrilateral":
-            # Handle quadrilateral short names like RTCF and RTCE.
-            element = element.reconstruct(cell=quadrilateral_tpc)
-        elif element.cell.cellname() == "hexahedron":
-            # Handle hexahedron short names like NCF and NCE.
-            element = element.reconstruct(cell=hexahedron_tpc)
-        else:
+        # Handle quadrilateral short names like RTCE/F and NCE/F.
+        try:
+            tpc = supported_tensor_product_cells[element.cell.cellname()]
+        except ValueError:
             raise ValueError("%s is supported, but handled incorrectly" %
                              element.family())
+        element = element.reconstruct(cell=tpc)
         finat_elem, deps = _create_element(element, **kwargs)
         return finat.FlattenedDimensions(finat_elem), deps
 
@@ -310,6 +312,11 @@ def convert_restrictedelement(element, **kwargs):
 
 hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval)
 quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval)
+supported_tensor_product_cells = {
+    "quadrilateral": quadrilateral_tpc,
+    "hexahedron": hexahedron_tpc,
+}
+
 _cache = weakref.WeakKeyDictionary()
 
 
diff --git a/finat/ufl/finiteelement.py b/finat/ufl/finiteelement.py
index 4cf34c90d..bec2854f3 100644
--- a/finat/ufl/finiteelement.py
+++ b/finat/ufl/finiteelement.py
@@ -94,6 +94,21 @@ def __new__(cls,
                 return EnrichedElement(HCurl(TensorProductElement(Qc_elt, Id_elt, cell=cell)),
                                        HCurl(TensorProductElement(Qd_elt, Ic_elt, cell=cell)))
 
+            elif family == "HDiv Trace":
+                cell_h, cell_v = cell.sub_cells()
+                cell_h, cell_v = cell.sub_cells()
+
+                hdegree = 0 if cell_h.cellname() == "interval" else degree
+                vdegree = 0 if cell_v.cellname() == "interval" else degree
+                tr_h = FiniteElement("HDiv Trace", cell_h, hdegree, variant=variant)
+                tr_v = FiniteElement("HDiv Trace", cell_v, vdegree, variant=variant)
+
+                dg_h = FiniteElement("DG", cell_h, degree, variant=variant)
+                dg_v = FiniteElement("DG", cell_v, degree, variant=variant)
+
+                return EnrichedElement(TensorProductElement(tr_h, dg_v, cell=cell),
+                                       TensorProductElement(dg_h, tr_v, cell=cell))
+
             elif family == "Q":
                 return TensorProductElement(*[FiniteElement("CG", c, degree, variant=variant)
                                               for c in cell.sub_cells()],
diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py
index 745518a75..d2586f57e 100644
--- a/test/FIAT/unit/test_fiat.py
+++ b/test/FIAT/unit/test_fiat.py
@@ -34,7 +34,7 @@
 from FIAT.nedelec import Nedelec                                # noqa: F401
 from FIAT.nedelec_second_kind import NedelecSecondKind          # noqa: F401
 from FIAT.regge import Regge                                    # noqa: F401
-from FIAT.hdiv_trace import HDivTrace, map_to_reference_facet   # noqa: F401
+from FIAT.hdiv_trace import HDivTrace                           # 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
@@ -400,7 +400,6 @@ def __init__(self, a, b):
     xfail_impl("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))"),
     xfail_impl("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"),
     xfail_impl("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"),
-    xfail_impl("HDivTrace(T, 1)"),
     xfail_impl("EnrichedElement("
                "Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
                "Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
@@ -409,15 +408,16 @@ def __init__(self, a, b):
                "Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
                "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
                ")"),
+    "HDivTrace(I, 0)",
+    "HDivTrace(T, 0)",
+    "HDivTrace(T, 1)",
+    "HDivTrace(T, 2)",
+    "HDivTrace(T, 3)",
+    "HDivTrace(S, 0)",
+    "HDivTrace(S, 1)",
+    "HDivTrace(S, 2)",
+    "HDivTrace(S, 3)",
     # Following elements are checked using tabulate
-    xfail_impl("HDivTrace(T, 0)"),
-    xfail_impl("HDivTrace(T, 1)"),
-    xfail_impl("HDivTrace(T, 2)"),
-    xfail_impl("HDivTrace(T, 3)"),
-    xfail_impl("HDivTrace(S, 0)"),
-    xfail_impl("HDivTrace(S, 1)"),
-    xfail_impl("HDivTrace(S, 2)"),
-    xfail_impl("HDivTrace(S, 3)"),
     xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"),
     xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"),
     xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"),
@@ -593,7 +593,6 @@ def test_facet_nodality_tabulate(element):
     facet_dim = sorted(entity_dofs.keys())[-2]
     facet_dofs = entity_dofs[facet_dim]
     dofs = element.dual_basis()
-    vertices = element.ref_el.vertices
 
     for (facet, indices) in facet_dofs.items():
         for i in indices:
@@ -602,16 +601,13 @@ def test_facet_nodality_tabulate(element):
             (coords, weights), = node.get_point_dict().items()
             assert weights == [(1.0, ())]
 
-            # Map dof coordinates to reference element due to
-            # HdivTrace interface peculiarity
-            ref_coords, = map_to_reference_facet((coords,), vertices, facet)
-            nodes_coords.append((facet, ref_coords))
+            nodes_coords.append(coords)
 
     # Check nodality
-    for j, (facet, x) in enumerate(nodes_coords):
-        basis, = element.tabulate(0, (x,), entity=(facet_dim, facet)).values()
-        for i in range(len(basis)):
-            assert np.isclose(basis[i], 1.0 if i == j else 0.0)
+    expected = np.eye(len(dofs))
+    for j, x in enumerate(nodes_coords):
+        basis, = element.tabulate(0, (x,)).values()
+        assert np.allclose(basis.T, expected[j])
 
 
 @pytest.mark.parametrize('element', [
diff --git a/test/FIAT/unit/test_hdivtrace.py b/test/FIAT/unit/test_hdivtrace.py
index ec405a190..9c5be7505 100644
--- a/test/FIAT/unit/test_hdivtrace.py
+++ b/test/FIAT/unit/test_hdivtrace.py
@@ -36,12 +36,13 @@ def test_basis_values(dim, degree, variant):
     (2) The entity pair (dim, id) is provided, and the trace element
         tabulates accordingly using the new tabulate API.
     """
-    from FIAT import ufc_simplex, HDivTrace, make_quadrature
+    from FIAT import ufc_simplex, HDivTrace, make_quadrature, DiscontinuousLagrange
 
     ref_el = ufc_simplex(dim)
-    quadrule = make_quadrature(ufc_simplex(dim - 1), degree + 1)
+    facet = ref_el.get_facet_element()
+    quadrule = make_quadrature(facet, degree + 1)
     fiat_element = HDivTrace(ref_el, degree, variant=variant)
-    facet_element = fiat_element.dg_elements[dim - 1]
+    facet_element = DiscontinuousLagrange(facet, degree, variant=variant)
     nf = facet_element.space_dimension()
 
     for facet_id in range(dim + 1):
@@ -70,20 +71,22 @@ def test_basis_values(dim, degree, variant):
 
 
 @pytest.mark.parametrize("degree", range(4))
-def test_quad_trace(degree):
+@pytest.mark.parametrize("variant", ("spectral", "integral"))
+def test_quad_trace(degree, variant):
     """Test the trace element defined on a quadrilateral cell"""
-    from FIAT import ufc_simplex, HDivTrace, make_quadrature
-    from FIAT.reference_element import TensorProductCell
+    from FIAT import ufc_cell, HDivTrace, make_quadrature, DiscontinuousLagrange
+
+    ref_el = ufc_cell("quadrilateral")
+    fiat_element = HDivTrace(ref_el, degree, variant=variant)
 
-    tpc = TensorProductCell(ufc_simplex(1), ufc_simplex(1))
-    fiat_element = HDivTrace(tpc, (degree, degree))
-    facet_elements = fiat_element.dg_elements
-    quadrule = make_quadrature(ufc_simplex(1), degree + 1)
+    facet = ref_el.get_facet_element()
+    facet_element = DiscontinuousLagrange(facet, degree, variant=variant)
+    quadrule = make_quadrature(facet, degree + 1)
 
-    for i, entity in enumerate([((0, 1), 0), ((0, 1), 1),
-                                ((1, 0), 0), ((1, 0), 1)]):
-        entity_dim, _ = entity
-        element = facet_elements[entity_dim]
+    entity_dim = facet.get_dimension()
+    for i, entity_id in enumerate(ref_el.topology[entity_dim]):
+        entity = (entity_dim, entity_id)
+        element = facet_element
         nf = element.space_dimension()
 
         tab = fiat_element.tabulate(0, quadrule.pts,

From 2d802c90bdf88842a87303c029205d799b53bb5b Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Fri, 4 Jul 2025 11:56:18 +0100
Subject: [PATCH 07/11] LagrangeExpansionSet on 1D manifolds

---
 FIAT/barycentric_interpolation.py | 38 ++++++++++++++++++++++++-------
 FIAT/discontinuous_lagrange.py    |  2 +-
 FIAT/expansions.py                | 16 ++++++-------
 finat/element_factory.py          |  2 +-
 4 files changed, 40 insertions(+), 18 deletions(-)

diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py
index 25862bf2f..4d19457a5 100644
--- a/FIAT/barycentric_interpolation.py
+++ b/FIAT/barycentric_interpolation.py
@@ -30,7 +30,7 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0):
         sp_simplify = numpy.vectorize(simplify)
     else:
         sp_simplify = lambda x: x
-    phi = numpy.add.outer(-nodes, pts.flatten())
+    phi = numpy.add.outer(-nodes.flatten(), pts.flatten())
     with numpy.errstate(divide='ignore', invalid='ignore'):
         numpy.reciprocal(phi, out=phi)
         numpy.multiply(phi, wts[:, None], out=phi)
@@ -49,7 +49,7 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0):
 def make_dmat(x):
     """Returns Lagrange differentiation matrix and barycentric weights
     associated with x[j]."""
-    dmat = numpy.add.outer(-x, x)
+    dmat = numpy.add.outer(-x.flatten(), x.flatten())
     numpy.fill_diagonal(dmat, 1.0)
     wts = numpy.prod(dmat, axis=0)
     numpy.reciprocal(wts, out=wts)
@@ -59,22 +59,38 @@ def make_dmat(x):
 
 
 class LagrangeLineExpansionSet(expansions.LineExpansionSet):
-    """Lagrange polynomial expansion set for given points the line."""
+    """Lagrange polynomial expansion set for given points on the line."""
     def __init__(self, ref_el, pts):
+        if ref_el.get_spatial_dimension() != 1:
+            raise Exception("Must have a line")
+        pts = numpy.asarray(pts)
         self.points = pts
-        self.x = numpy.array(pts, dtype="d").flatten()
+
         self.cell_node_map = expansions.compute_cell_point_map(ref_el, pts, unique=False)
         self.dmats = [None for _ in self.cell_node_map]
         self.weights = [None for _ in self.cell_node_map]
         self.nodes = [None for _ in self.cell_node_map]
+        self.affine_mappings = {}
         for cell, ibfs in self.cell_node_map.items():
-            self.nodes[cell] = self.x[ibfs]
-            self.dmats[cell], self.weights[cell] = make_dmat(self.nodes[cell])
+            x = pts[ibfs]
+            if ref_el.is_trace():
+                verts = ref_el.get_vertices_of_subcomplex(ref_el.topology[1][cell])
+                A = numpy.diff(verts, axis=0)[0]/2
+                A /= numpy.linalg.norm(A)
+                b = -numpy.dot(numpy.sum(verts, axis=0)/2, A.T)
+                self.affine_mappings[cell] = (A, b)
+                x = numpy.add(numpy.dot(x, A.T), b)
+            self.nodes[cell] = x
+            self.dmats[cell], self.weights[cell] = make_dmat(x)
 
         self.degree = max(len(wts) for wts in self.weights)-1
         self.recurrence_order = self.degree + 1
-        super().__init__(ref_el)
-        self.continuity = None if len(self.x) == sum(len(xk) for xk in self.nodes) else "C0"
+        self.ref_el = ref_el
+        self.variant = None
+        self.scale = 1
+        self.continuity = None if len(pts) == sum(len(xk) for xk in self.nodes) else "C0"
+        self._dmats_cache = {}
+        self._cell_node_map_cache = {}
 
     def get_num_members(self, n):
         return len(self.points)
@@ -89,6 +105,12 @@ def get_dmats(self, degree, cell=0):
         return [self.dmats[cell].T]
 
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
+        try:
+            A, b = self.affine_mappings[cell]
+            ref_pts = numpy.add(numpy.dot(pts.reshape(-1, A.shape[-1]), A.T), b)
+            pts = ref_pts.reshape(*pts.shape[:-1], -1)
+        except KeyError:
+            pass
         return barycentric_interpolation(self.nodes[cell], self.weights[cell], self.dmats[cell], pts, order=order)
 
 
diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py
index db5c06b83..f806c9278 100644
--- a/FIAT/discontinuous_lagrange.py
+++ b/FIAT/discontinuous_lagrange.py
@@ -233,7 +233,7 @@ def __init__(self, ref_el, degree, variant="equispaced"):
             dual = BrokenLagrangeDualSet(ref_el, degree, point_variant=point_variant)
         else:
             dual = DiscontinuousLagrangeDualSet(ref_el, degree, point_variant=point_variant)
-        if ref_el.shape == LINE and not ref_el.is_trace():
+        if ref_el.shape == LINE:
             # In 1D we can use the primal basis as the expansion set,
             # avoiding any round-off coming from a basis transformation
             points = get_lagrange_points(dual)
diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index 959490552..100163e29 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -371,7 +371,7 @@ def _tabulate(self, n, pts, order=0):
             if tdim == 0 and len(phis) == 0:
                 # Hack for TensorProduct HDivTrace: do not raise TraceError on the interval
                 for cell in parent.topology[gdim]:
-                    phis[cell] = {(0,)*gdim: numpy.zeros((1, 1))}
+                    phis[cell] = {(0,)*gdim: numpy.zeros(())}
             elif sum(len(cell_point_map[cell]) for cell in cell_point_map) < len(pts):
                 # Raise TraceError when interior points fail to be binned on facets
                 for cell in parent.topology[gdim]:
@@ -626,14 +626,14 @@ def get_affine_mapping(xs, ys):
     """
     X = numpy.asarray(xs)
     Y = numpy.asarray(ys)
-    A = X[1:] - X[:1]
-    B = Y[1:] - Y[:1]
-    if A.shape[0] == A.shape[1]:
-        C = numpy.linalg.solve(A, B)
+    DX = X[1:] - X[:1]
+    DY = Y[1:] - Y[:1]
+    if DX.shape[0] == DX.shape[1]:
+        AT = numpy.linalg.solve(DX, DY)
     else:
-        C, *_ = numpy.linalg.lstsq(A, B)
-    b = Y[0] - numpy.dot(X[0], C)
-    return C.T, b
+        AT, *_ = numpy.linalg.lstsq(DX, DY)
+    b = Y[0] - numpy.dot(X[0], AT)
+    return AT.T, b
 
 
 def polynomial_dimension(ref_el, n, continuity=None):
diff --git a/finat/element_factory.py b/finat/element_factory.py
index fcca67d3e..d0d457d49 100644
--- a/finat/element_factory.py
+++ b/finat/element_factory.py
@@ -169,7 +169,7 @@ def convert_finiteelement(element, **kwargs):
         # Handle quadrilateral short names like RTCE/F and NCE/F.
         try:
             tpc = supported_tensor_product_cells[element.cell.cellname()]
-        except ValueError:
+        except KeyError:
             raise ValueError("%s is supported, but handled incorrectly" %
                              element.family())
         element = element.reconstruct(cell=tpc)

From c0597596aeb3aad44e7153b7cf06929d72bb3b63 Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Fri, 4 Jul 2025 17:55:56 +0100
Subject: [PATCH 08/11] finat.ufl.HDivTraceElement, full support for tensor
 product elements

---
 FIAT/barycentric_interpolation.py |  2 +-
 FIAT/expansions.py                | 57 +++++++++++--------
 FIAT/hdiv_trace.py                |  7 +--
 FIAT/reference_element.py         | 11 ++--
 finat/__init__.py                 |  2 +-
 finat/element_factory.py          |  6 ++
 finat/fiat_elements.py            |  3 +-
 finat/hdivcurl.py                 | 60 +++++++++++++++++++-
 finat/ufl/__init__.py             |  2 +-
 finat/ufl/finiteelement.py        | 22 +++++---
 finat/ufl/hdivcurl.py             | 93 ++++++++++---------------------
 11 files changed, 151 insertions(+), 114 deletions(-)

diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py
index 4d19457a5..8b1fba4a4 100644
--- a/FIAT/barycentric_interpolation.py
+++ b/FIAT/barycentric_interpolation.py
@@ -75,7 +75,7 @@ def __init__(self, ref_el, pts):
             x = pts[ibfs]
             if ref_el.is_trace():
                 verts = ref_el.get_vertices_of_subcomplex(ref_el.topology[1][cell])
-                A = numpy.diff(verts, axis=0)[0]/2
+                A, = numpy.diff(verts, axis=0)
                 A /= numpy.linalg.norm(A)
                 b = -numpy.dot(numpy.sum(verts, axis=0)/2, A.T)
                 self.affine_mappings[cell] = (A, b)
diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index 100163e29..c0780b324 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -345,8 +345,6 @@ def distance(alpha, beta):
 
     def _tabulate(self, n, pts, order=0):
         """A version of tabulate() that also works for a single point."""
-        from FIAT.hdiv_trace import TraceError
-        from FIAT.polynomial_set import mis
         pts = numpy.asarray(pts)
         unique = self.continuity is not None and order == 0
         cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=unique)
@@ -357,26 +355,7 @@ def _tabulate(self, n, pts, order=0):
             return phis[0]
 
         if self.ref_el.is_trace():
-            parent = self.ref_el.get_parent()
-            tdim = self.ref_el.get_spatial_dimension()
-            gdim = parent.get_spatial_dimension()
-            for cell in phis:
-                # Promote facet keys to cell keys
-                phi = phis[cell][(0,) * tdim]
-                # Raise TraceError on gradient tabulations
-                msg = "Gradients on trace elements are not well-defined."
-                phis[cell] = {alpha: phi if sum(alpha) == 0 else TraceError(msg)
-                              for i in range(order+1)
-                              for alpha in mis(gdim, i)}
-            if tdim == 0 and len(phis) == 0:
-                # Hack for TensorProduct HDivTrace: do not raise TraceError on the interval
-                for cell in parent.topology[gdim]:
-                    phis[cell] = {(0,)*gdim: numpy.zeros(())}
-            elif sum(len(cell_point_map[cell]) for cell in cell_point_map) < len(pts):
-                # Raise TraceError when interior points fail to be binned on facets
-                for cell in parent.topology[gdim]:
-                    msg = "The HDivTrace element can only be tabulated on facets."
-                    phis[cell] = {(0,)*gdim: TraceError(msg)}
+            phis = trace_tabulation(self.ref_el, cell_point_map, order, pts, phis)
 
         if pts.dtype == object:
             # If binning is undefined, scale by the characteristic function of each subcell
@@ -402,7 +381,7 @@ def _tabulate(self, n, pts, order=0):
         result = {}
         base_phi = tuple(phis.values())[0]
         for alpha in base_phi:
-            if isinstance(base_phi[alpha], TraceError):
+            if isinstance(base_phi[alpha], Exception):
                 result[alpha] = base_phi[alpha]
                 continue
             dtype = base_phi[alpha].dtype
@@ -723,7 +702,10 @@ def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12):
         return {cell: Ellipsis for cell in sorted(top[sd])}
 
     # The distance to the nearest cell is equal to the distance to the parent cell
-    best = ref_el.get_parent().distance_to_point_l1(pts, rescale=True)
+    parent = ref_el
+    while parent.get_parent() is not None:
+        parent = parent.get_parent()
+    best = parent.distance_to_point_l1(pts, rescale=True)
     tol = best + tol
 
     cell_point_map = {}
@@ -779,3 +761,30 @@ def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12):
         mult = sum(masks)
         masks = [m / mult for m in masks]
     return masks
+
+
+def trace_tabulation(ref_el, cell_point_map, order, pts, phis):
+    """Lift trace tabulations into the cells and raise TraceError on invalid tabulations."""
+    from FIAT.polynomial_set import mis
+    from FIAT.hdiv_trace import TraceError
+    parent = ref_el.get_parent()
+    tdim = ref_el.get_spatial_dimension()
+    gdim = parent.get_spatial_dimension()
+    facet_key = (0, ) * tdim
+    cell_key = (0, ) * gdim
+
+    for cell in phis:
+        # Lift facet keys to cell keys
+        phi = phis[cell][facet_key]
+        # Raise TraceError on gradient tabulations
+        msg = "Gradients on trace elements are not well-defined."
+        phis[cell] = {alpha: phi if sum(alpha) == 0 else TraceError(msg)
+                      for i in range(order+1)
+                      for alpha in mis(gdim, i)}
+
+    if sum(len(cell_point_map[cell]) for cell in cell_point_map) < len(pts):
+        # Raise TraceError when interior points fail to be binned on facets
+        for cell in parent.topology[gdim]:
+            msg = "The HDivTrace element can only be tabulated on facets."
+            phis[cell] = {cell_key: TraceError(msg)}
+    return phis
diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py
index 2ea70a4ee..3c26d0b47 100644
--- a/FIAT/hdiv_trace.py
+++ b/FIAT/hdiv_trace.py
@@ -100,11 +100,8 @@ class HDivTrace(DiscontinuousLagrange):
     def __new__(cls, ref_el, degree, variant="equispaced_interior"):
         """Constructor for the HDivTrace element.
 
-        :arg ref_el: A reference element, which may be a tensor product
-                     cell.
-        :arg degree: The degree of approximation. If on a tensor product
-                     cell, then provide a tuple of degrees if you want
-                     varying degrees.
+        :arg ref_el: A reference element.
+        :arg degree: The degree of approximation.
         :arg variant: The point distribution variant passed on to recursivenodes.
         """
         facets = TraceSimplicialComplex(ref_el)
diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index 8a2aed15e..9744e0bd0 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -618,7 +618,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
         if len(points) == 0:
             return points
         if entity is None:
-            entity = (self.get_spatial_dimension(), 0)
+            entity = (self.get_dimension(), 0)
         entity_dim, entity_id = entity
         sd = len(self.vertices[0])
 
@@ -628,10 +628,9 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
         top = self.get_topology()
         subcomplex = top[entity_dim][entity_id]
         if entity_dim != sd:
-            parent = self.get_parent() if self.is_trace() else self
-            cell_id = parent.connectivity[(entity_dim, sd)][0][0]
-            top = parent.get_topology()
-            subcell = top[sd][cell_id]
+            parent = self.get_parent_complex() if self.is_trace() else self
+            cell_id = min(parent.connectivity[(entity_dim, sd)][entity_id])
+            subcell = parent.topology[sd][cell_id]
             while len(subcell) > sd + 1:
                 # construct a simplex if we have a hypercube
                 k = max(set(subcell) - set(subcomplex))
@@ -648,7 +647,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
         A, b = make_affine_mapping(cell_verts, ref_verts)
         A, b = A[indices], b[indices]
         if rescale:
-            # rescale barycentric coordinates by the height wrt. to the facet
+            # rescale barycentric coordinates by the height w.r.t. the facet
             h = 1 / numpy.linalg.norm(A, axis=1)
             b *= h
             A *= h[:, None]
diff --git a/finat/__init__.py b/finat/__init__.py
index 5390fa779..dde090cb0 100644
--- a/finat/__init__.py
+++ b/finat/__init__.py
@@ -32,7 +32,7 @@
 from .cube import FlattenedDimensions                              # noqa: F401
 from .discontinuous import DiscontinuousElement                    # noqa: F401
 from .enriched import EnrichedElement                              # noqa: F401
-from .hdivcurl import HCurlElement, HDivElement                    # noqa: F401
+from .hdivcurl import HCurlElement, HDivElement, HDivTraceElement  # noqa: F401
 from .mixed import MixedElement                                    # noqa: F401
 from .nodal_enriched import NodalEnrichedElement                   # noqa: F401
 from .quadrature_element import QuadratureElement, make_quadrature_element  # noqa: F401
diff --git a/finat/element_factory.py b/finat/element_factory.py
index d0d457d49..bb665849d 100644
--- a/finat/element_factory.py
+++ b/finat/element_factory.py
@@ -287,6 +287,12 @@ def convert_tensorproductelement(element, **kwargs):
     return finat.TensorProductElement(elements), deps
 
 
+@convert.register(finat.ufl.HDivTraceElement)
+def convert_hdivtraceelement(element, **kwargs):
+    finat_elem, deps = _create_element(element._element, **kwargs)
+    return finat.HDivTraceElement(finat_elem), deps
+
+
 @convert.register(finat.ufl.HDivElement)
 def convert_hdivelement(element, **kwargs):
     finat_elem, deps = _create_element(element._element, **kwargs)
diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py
index 3677ef965..b94cc1a7e 100644
--- a/finat/fiat_elements.py
+++ b/finat/fiat_elements.py
@@ -112,7 +112,8 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         index_shape = (self._element.space_dimension(),)
         for alpha, fiat_table in fiat_result.items():
             if isinstance(fiat_table, Exception):
-                result[alpha] = gem.Failure(self.index_shape + self.value_shape, fiat_table)
+                shape = ps.points.shape[:-1] + self.index_shape + self.value_shape
+                result[alpha] = gem.partial_indexed(gem.Failure(shape, fiat_table), ps.indices)
                 continue
 
             derivative = sum(alpha)
diff --git a/finat/hdivcurl.py b/finat/hdivcurl.py
index 3ca697744..c0bdb4454 100644
--- a/finat/hdivcurl.py
+++ b/finat/hdivcurl.py
@@ -2,6 +2,7 @@
 from FIAT.reference_element import LINE
 
 import gem
+import numpy
 from gem.utils import cached_property
 from finat.finiteelementbase import FiniteElementBase
 from finat.tensor_product import TensorProductElement
@@ -58,7 +59,7 @@ def index_shape(self):
     def value_shape(self):
         return (self.cell.get_spatial_dimension(),)
 
-    def _transform_evaluation(self, core_eval):
+    def _transform_evaluation(self, core_eval, entity=None):
         beta = self.get_indices()
         zeta = self.get_value_indices()
 
@@ -72,11 +73,11 @@ def promote(table):
 
     def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
         core_eval = self.wrappee.basis_evaluation(order, ps, entity)
-        return self._transform_evaluation(core_eval)
+        return self._transform_evaluation(core_eval, entity)
 
     def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None):
         core_eval = self.wrappee.point_evaluation(order, refcoords, entity)
-        return self._transform_evaluation(core_eval)
+        return self._transform_evaluation(core_eval, entity)
 
     @property
     def dual_basis(self):
@@ -92,6 +93,59 @@ def dual_basis(self):
         return gem.ComponentTensor(Q[zeta], beta + zeta), x
 
 
+class HDivTraceElement(WrapperElementBase):
+    """HDiv Trace wrapper element for tensor product elements."""
+
+    def __init__(self, wrappee):
+        assert isinstance(wrappee, TensorProductElement)
+        if any(fe.formdegree is None for fe in wrappee.factors):
+            raise ValueError("Form degree of subelement is None, cannot HDiv Trace!")
+
+        formdegree = sum(fe.formdegree for fe in wrappee.factors)
+        if formdegree != wrappee.cell.get_spatial_dimension() - 1:
+            raise ValueError("HDiv Trace requires (n-1)-form element!")
+
+        self.support_dim = tuple(fe.formdegree for fe in wrappee.factors)
+
+        super().__init__(wrappee, None)
+
+    def _transform_evaluation(self, core_eval, entity):
+        if entity is None:
+            entity = (self.cell.get_dimension(), 0)
+        entity_dim, entity_id = entity
+
+        if entity_dim == self.support_dim or sum(entity_dim) != sum(self.support_dim):
+            return core_eval
+
+        def zero_failure(expr):
+            if isinstance(expr, gem.Failure):
+                return gem.Literal(numpy.zeros(expr.shape, expr.dtype))
+            return expr.reconstruct(*map(zero_failure, expr.children))
+
+        # Create a zero tabulation with the same tensor-product structure
+        return {alpha: zero_failure(table) for alpha, table in core_eval.items()}
+
+    @property
+    def formdegree(self):
+        return self.wrappee.formdegree
+
+    @property
+    def value_shape(self):
+        return self.wrappee.value_shape
+
+    @cached_property
+    def fiat_equivalent(self):
+        return self.wrappee.fiat_equivalent
+
+    @property
+    def mapping(self):
+        return self.wrappee.mapping
+
+    @property
+    def dual_basis(self):
+        return self.wrappee.dual_basis
+
+
 class HDivElement(WrapperElementBase):
     """H(div) wrapper element for tensor product elements."""
 
diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py
index 21a7d13d1..f36ca39db 100644
--- a/finat/ufl/__init__.py
+++ b/finat/ufl/__init__.py
@@ -15,7 +15,7 @@
 from finat.ufl.enrichedelement import EnrichedElement, NodalEnrichedElement  # noqa: F401
 from finat.ufl.finiteelement import FiniteElement  # noqa: F401
 from finat.ufl.finiteelementbase import FiniteElementBase  # noqa: F401
-from finat.ufl.hdivcurl import HCurlElement, HDivElement, WithMapping, HDiv, HCurl  # noqa: F401
+from finat.ufl.hdivcurl import HCurlElement, HDivElement, HDivTraceElement, WithMapping, HDiv, HCurl  # noqa: F401
 from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement  # noqa: F401
 from finat.ufl.restrictedelement import RestrictedElement  # noqa: F401
 from finat.ufl.tensorproductelement import TensorProductElement  # noqa: F401
diff --git a/finat/ufl/finiteelement.py b/finat/ufl/finiteelement.py
index bec2854f3..8a2103665 100644
--- a/finat/ufl/finiteelement.py
+++ b/finat/ufl/finiteelement.py
@@ -39,6 +39,7 @@ def __new__(cls,
             from finat.ufl.enrichedelement import EnrichedElement
             from finat.ufl.hdivcurl import HCurlElement as HCurl
             from finat.ufl.hdivcurl import HDivElement as HDiv
+            from finat.ufl.hdivcurl import HDivTraceElement as HDivTrace
             from finat.ufl.tensorproductelement import TensorProductElement
 
             family, short_name, degree, reference_value_shape, sobolev_space, mapping, embedded_degree = \
@@ -96,18 +97,21 @@ def __new__(cls,
 
             elif family == "HDiv Trace":
                 cell_h, cell_v = cell.sub_cells()
-                cell_h, cell_v = cell.sub_cells()
+                if not isinstance(degree, tuple):
+                    degree = (degree, degree)
+                hdegree, vdegree = degree
+
+                dg_family = lambda cell: "DG" if cell.is_simplex() else "DQ"
+                is_interval = lambda cell: cell.cellname() == "interval"
 
-                hdegree = 0 if cell_h.cellname() == "interval" else degree
-                vdegree = 0 if cell_v.cellname() == "interval" else degree
-                tr_h = FiniteElement("HDiv Trace", cell_h, hdegree, variant=variant)
-                tr_v = FiniteElement("HDiv Trace", cell_v, vdegree, variant=variant)
+                dg_h = FiniteElement(dg_family(cell_h), cell_h, hdegree, variant=variant)
+                tr_h = FiniteElement("HDiv Trace", cell_h, 0 if is_interval(cell_h) else hdegree, variant=variant)
 
-                dg_h = FiniteElement("DG", cell_h, degree, variant=variant)
-                dg_v = FiniteElement("DG", cell_v, degree, variant=variant)
+                dg_v = FiniteElement(dg_family(cell_v), cell_v, vdegree, variant=variant)
+                tr_v = FiniteElement("HDiv Trace", cell_v, 0 if is_interval(cell_v) else vdegree, variant=variant)
 
-                return EnrichedElement(TensorProductElement(tr_h, dg_v, cell=cell),
-                                       TensorProductElement(dg_h, tr_v, cell=cell))
+                return EnrichedElement(HDivTrace(TensorProductElement(tr_h, dg_v, cell=cell)),
+                                       HDivTrace(TensorProductElement(dg_h, tr_v, cell=cell)))
 
             elif family == "Q":
                 return TensorProductElement(*[FiniteElement("CG", c, degree, variant=variant)
diff --git a/finat/ufl/hdivcurl.py b/finat/ufl/hdivcurl.py
index ed188038e..0a80e397d 100644
--- a/finat/ufl/hdivcurl.py
+++ b/finat/ufl/hdivcurl.py
@@ -35,40 +35,38 @@ def __call__(self, element):
 HDiv = CallableSobolevSpace(HDivSobolevSpace.name, HDivSobolevSpace.parents)
 
 
-class HDivElement(FiniteElementBase):
-    """A div-conforming version of an outer product element, assuming this makes mathematical sense."""
+class WrapperElementBase(FiniteElementBase):
+    """A modified version of a tensor product element."""
     __slots__ = ("_element", )
 
-    def __init__(self, element):
-        """Doc."""
+    def __init__(self, element, reference_value_shape, sobolev_space, mapping):
         self._element = element
+        self._sobolev_space = sobolev_space
+        self._mapping = mapping
 
-        family = "TensorProductElement"
+        family = element.family()
         cell = element.cell
         degree = element.degree()
         quad_scheme = element.quadrature_scheme()
-        reference_value_shape = (element.cell.topological_dimension(),)
 
-        # Skipping TensorProductElement constructor! Bad code smell, refactor to avoid this non-inheritance somehow.
-        FiniteElementBase.__init__(self, family, cell, degree,
-                                   quad_scheme, reference_value_shape)
+        super().__init__(family, cell, degree, quad_scheme, reference_value_shape)
 
     def __repr__(self):
         """Doc."""
-        return f"HDivElement({repr(self._element)})"
+        return f"{type(self).__name__}({repr(self._element)})"
 
     def mapping(self):
         """Doc."""
-        return "contravariant Piola"
+        return self._mapping
 
     @property
     def sobolev_space(self):
         """Return the underlying Sobolev space."""
-        return HDivSobolevSpace
+        return self._sobolev_space
 
     def reconstruct(self, **kwargs):
         """Doc."""
-        return HDivElement(self._element.reconstruct(**kwargs))
+        return type(self)(self._element.reconstruct(**kwargs))
 
     def variant(self):
         """Doc."""
@@ -76,11 +74,11 @@ def variant(self):
 
     def __str__(self):
         """Doc."""
-        return f"HDivElement({repr(self._element)})"
+        return f"{type(self).__name__}({repr(self._element)})"
 
     def shortstr(self):
         """Format as string for pretty printing."""
-        return f"HDivElement({self._element.shortstr()})"
+        return f"HDivTraceElement({self._element.shortstr()})"
 
     @property
     def embedded_subdegree(self):
@@ -93,63 +91,32 @@ def embedded_superdegree(self):
         return self._element.embedded_superdegree
 
 
-class HCurlElement(FiniteElementBase):
-    """A curl-conforming version of an outer product element, assuming this makes mathematical sense."""
-    __slots__ = ("_element",)
+class HDivTraceElement(WrapperElementBase):
+    """A HDiv Trace version of a tensor product element, assuming this makes mathematical sense."""
+    __slots__ = ("_element", )
 
     def __init__(self, element):
-        """Doc."""
-        self._element = element
+        reference_value_shape = ()
+        super().__init__(element, reference_value_shape, L2, "identity")
 
-        family = "TensorProductElement"
-        cell = element.cell
-        degree = element.degree()
-        quad_scheme = element.quadrature_scheme()
-        cell = element.cell
-        reference_value_shape = (cell.topological_dimension(),)  # TODO: Is this right?
-        # Skipping TensorProductElement constructor! Bad code smell,
-        # refactor to avoid this non-inheritance somehow.
-        FiniteElementBase.__init__(self, family, cell, degree, quad_scheme,
-                                   reference_value_shape)
-
-    def __repr__(self):
-        """Doc."""
-        return f"HCurlElement({repr(self._element)})"
 
-    def mapping(self):
-        """Doc."""
-        return "covariant Piola"
+class HDivElement(WrapperElementBase):
+    """A div-conforming version of a tensor product element, assuming this makes mathematical sense."""
+    __slots__ = ("_element", )
 
-    @property
-    def sobolev_space(self):
-        """Return the underlying Sobolev space."""
-        return HCurlSobolevSpace
+    def __init__(self, element):
+        reference_value_shape = (element.cell.topological_dimension(),)
+        super().__init__(element, reference_value_shape, HDiv, "contravariant Piola")
 
-    def reconstruct(self, **kwargs):
-        """Doc."""
-        return HCurlElement(self._element.reconstruct(**kwargs))
 
-    def variant(self):
-        """Doc."""
-        return self._element.variant()
+class HCurlElement(WrapperElementBase):
+    """A curl-conforming version of a tensor product element, assuming this makes mathematical sense."""
+    __slots__ = ("_element",)
 
-    def __str__(self):
+    def __init__(self, element):
         """Doc."""
-        return f"HCurlElement({repr(self._element)})"
-
-    def shortstr(self):
-        """Format as string for pretty printing."""
-        return f"HCurlElement({self._element.shortstr()})"
-
-    @property
-    def embedded_subdegree(self):
-        """Return embedded subdegree."""
-        return self._element.embedded_subdegree
-
-    @property
-    def embedded_superdegree(self):
-        """Return embedded superdegree."""
-        return self._element.embedded_superdegree
+        reference_value_shape = (element.cell.topological_dimension(),)
+        super().__init__(element, reference_value_shape, HCurl, "covariant Piola")
 
 
 class WithMapping(FiniteElementBase):

From 2769a35dfb9537d174a51319a8153317cd3a71cd Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Mon, 7 Jul 2025 12:59:41 +0100
Subject: [PATCH 09/11] cleanup

---
 FIAT/reference_element.py | 37 +++++++++++++------------------------
 1 file changed, 13 insertions(+), 24 deletions(-)

diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index 9744e0bd0..553658acd 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -620,32 +620,21 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
         if entity is None:
             entity = (self.get_dimension(), 0)
         entity_dim, entity_id = entity
-        sd = len(self.vertices[0])
-
-        # get a subcell containing the entity and the restriction indices of the entity
-        outside = None
-        indices = slice(None)
-        top = self.get_topology()
-        subcomplex = top[entity_dim][entity_id]
-        if entity_dim != sd:
-            parent = self.get_parent_complex() if self.is_trace() else self
-            cell_id = min(parent.connectivity[(entity_dim, sd)][entity_id])
-            subcell = parent.topology[sd][cell_id]
-            while len(subcell) > sd + 1:
-                # construct a simplex if we have a hypercube
-                k = max(set(subcell) - set(subcomplex))
-                subcell = subcell[:k] + subcell[k+1:]
-
-            indices = [i for i, v in enumerate(subcell) if v in subcomplex]
-            subcomplex = subcell
+        subcomplex = self.topology[entity_dim][entity_id]
 
-            outside = [i for i in range(len(subcomplex)) if i not in indices]
-            indices.extend(outside)
+        tdim = sum(entity_dim) if isinstance(entity_dim, tuple) else entity_dim
+        gdim = len(self.vertices[0])
+        if tdim != gdim:
+            # Add extra vertices to construct a subcell containing the entity
+            parent = self.get_parent_complex() if self.is_trace() else self
+            cell_dim = parent.get_dimension()
+            cell_id = min(parent.connectivity[(entity_dim, cell_dim)][entity_id])
+            subcell = parent.topology[cell_dim][cell_id]
+            subcomplex += tuple(set(subcell) - set(subcomplex))[:gdim-tdim]
 
         cell_verts = self.get_vertices_of_subcomplex(subcomplex)
         ref_verts = numpy.eye(len(cell_verts))
         A, b = make_affine_mapping(cell_verts, ref_verts)
-        A, b = A[indices], b[indices]
         if rescale:
             # rescale barycentric coordinates by the height w.r.t. the facet
             h = 1 / numpy.linalg.norm(A, axis=1)
@@ -653,9 +642,9 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False):
             A *= h[:, None]
         bary = numpy.dot(points, A.T)
         numpy.add(bary, b, out=bary)
-
-        if outside is not None:
-            bary[:, -len(outside):] = -abs(bary[:, -len(outside):])
+        if tdim != gdim:
+            # the extra coordinates are always negated to mark points outside the facet
+            bary[:, tdim-gdim:] = -abs(bary[:, tdim-gdim:])
         return bary
 
     def compute_bubble(self, points, entity=None):

From 319354a4678b58a402e37ec6182e5174df88bc37 Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Tue, 8 Jul 2025 09:47:27 +0100
Subject: [PATCH 10/11] Add family to EnrichedElement

---
 finat/ufl/enrichedelement.py | 9 ++++++---
 finat/ufl/finiteelement.py   | 2 +-
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/finat/ufl/enrichedelement.py b/finat/ufl/enrichedelement.py
index a054885bb..bf686c694 100644
--- a/finat/ufl/enrichedelement.py
+++ b/finat/ufl/enrichedelement.py
@@ -17,7 +17,7 @@
 class EnrichedElementBase(FiniteElementBase):
     """The vector sum of several finite element spaces."""
 
-    def __init__(self, *elements):
+    def __init__(self, *elements, family=None):
         """Doc."""
         self._elements = elements
 
@@ -49,9 +49,11 @@ def __init__(self, *elements):
 
         # Get name of subclass: EnrichedElement or NodalEnrichedElement
         class_name = self.__class__.__name__
+        if family is None:
+            family = class_name
 
         # Initialize element data
-        FiniteElementBase.__init__(self, class_name, cell, degree,
+        FiniteElementBase.__init__(self, family, cell, degree,
                                    quad_scheme, reference_value_shape)
 
     def mapping(self):
@@ -87,7 +89,8 @@ def variant(self):
 
     def reconstruct(self, **kwargs):
         """Doc."""
-        return type(self)(*[e.reconstruct(**kwargs) for e in self._elements])
+        family = kwargs.pop("family", self.family())
+        return type(self)(*(e.reconstruct(**kwargs) for e in self._elements), family=family)
 
     @property
     def embedded_subdegree(self):
diff --git a/finat/ufl/finiteelement.py b/finat/ufl/finiteelement.py
index 8a2103665..0968f7e18 100644
--- a/finat/ufl/finiteelement.py
+++ b/finat/ufl/finiteelement.py
@@ -111,7 +111,7 @@ def __new__(cls,
                 tr_v = FiniteElement("HDiv Trace", cell_v, 0 if is_interval(cell_v) else vdegree, variant=variant)
 
                 return EnrichedElement(HDivTrace(TensorProductElement(tr_h, dg_v, cell=cell)),
-                                       HDivTrace(TensorProductElement(dg_h, tr_v, cell=cell)))
+                                       HDivTrace(TensorProductElement(dg_h, tr_v, cell=cell)), family=family)
 
             elif family == "Q":
                 return TensorProductElement(*[FiniteElement("CG", c, degree, variant=variant)

From 42bbd8063fcde21cf0b21de678e18d285333a565 Mon Sep 17 00:00:00 2001
From: Pablo Brubeck <brubeck@protonmail.com>
Date: Tue, 8 Jul 2025 13:57:26 +0100
Subject: [PATCH 11/11] Fixup Trace tabulation

---
 FIAT/expansions.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index c0780b324..fd1de72b8 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -770,12 +770,9 @@ def trace_tabulation(ref_el, cell_point_map, order, pts, phis):
     parent = ref_el.get_parent()
     tdim = ref_el.get_spatial_dimension()
     gdim = parent.get_spatial_dimension()
-    facet_key = (0, ) * tdim
-    cell_key = (0, ) * gdim
-
     for cell in phis:
         # Lift facet keys to cell keys
-        phi = phis[cell][facet_key]
+        phi = phis[cell][(0,) * tdim]
         # Raise TraceError on gradient tabulations
         msg = "Gradients on trace elements are not well-defined."
         phis[cell] = {alpha: phi if sum(alpha) == 0 else TraceError(msg)
@@ -786,5 +783,7 @@ def trace_tabulation(ref_el, cell_point_map, order, pts, phis):
         # Raise TraceError when interior points fail to be binned on facets
         for cell in parent.topology[gdim]:
             msg = "The HDivTrace element can only be tabulated on facets."
-            phis[cell] = {cell_key: TraceError(msg)}
+            phis[cell] = {alpha: TraceError(msg)
+                          for i in range(order+1)
+                          for alpha in mis(gdim, i)}
     return phis