diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 42984a86b..4e8ac8f81 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -16,7 +16,7 @@ jobs:
     strategy:
       matrix:
         os: [ubuntu-latest, macos-latest]
-        python-version: ['3.10', '3.11', '3.12', '3.13']
+        python-version: ['3.11', '3.12', '3.13', '3.14']
 
     steps:
       - uses: actions/checkout@v4
diff --git a/AGENTS.md b/AGENTS.md
new file mode 100644
index 000000000..af5952784
--- /dev/null
+++ b/AGENTS.md
@@ -0,0 +1,99 @@
+# AGENTS.md for FIAT
+
+This document outlines the guidelines and architectural context for AI agents assisting with the FIAT codebase, functioning as a core component within the broader Firedrake ecosystem. 
+
+---
+
+## AI Contribution Policy
+
+When assisting with contributions to FIAT and the Firedrake project, AI agents and their human counterparts must adhere to the following strict policies:
+
+* The use of AI tools must be explicitly declared alongside the specific tool used.
+* A human developer must lead the Pull Request.
+* The human contributor must understand every change made to the codebase.
+* Reviewer questions must be answered directly by the human, rather than acting as a relay to the AI.
+* Any generated code must be executed locally to verify that it functions correctly.
+* AI tools must not be used to resolve issues that are labeled as 'good first issue'.
+
+---
+
+## FIAT's Role in the Architecture
+
+FIAT operates within Firedrake's automated system for solving partial differential equations via the finite element method. Its specific architectural responsibilities include:
+
+* FIAT provides compile-time pre-tabulated basis functions.
+* These basis functions are utilized when the Two-Stage Form Compiler (TSFC) lowers Unified Form Language (UFL) into the GEM tensor language.
+* The resulting GEM expressions represent mathematical operations over quadrature points.
+* Mesh-topology bookkeeping routines within the ecosystem rely on FIAT ordering, such as the `create_cell_closure()` loop which builds a FIAT-ordered closure map necessary for subsequent code generation.
+
+---
+
+## Environment and Setup
+
+Bugs can exist within Firedrake or any of its component packages, explicitly including FIAT. To effectively develop and debug FIAT:
+
+* Developers should use editable installs for subpackages like FIAT so that source code edits take effect without requiring a full reinstallation.
+* The active branch or commit of each component must be verified before assuming a bug originates in the top-level Firedrake package.
+
+---
+
+## Core Coding Rules
+
+Agents modifying FIAT code must follow these fundamental development principles:
+
+* Bug fixes must target the underlying mathematical or architectural root cause.
+* Developers must avoid merely patching specific failing test cases or edge cases.
+* Code complexity should be minimized by favoring the mathematical generality of finite elements over complicated special-case logic.
+* Memorized API shapes must not be trusted.
+* APIs across the ecosystem evolve, meaning properties can become methods, arguments can be renamed, and signatures can be deprecated.
+* Agents must verify current API signatures by reading the installed source code instead of relying on outdated training data.
+* Code documentation and comments must explain the present, correct code.
+* Comments must not detail what a removed or incorrect approach previously did.
+
+---
+
+## Pattern Matching and Mathematical Reasoning
+
+When designing or debugging FIAT, FInAT, and GEM changes, use the existing codebase as a library of
+mathematical patterns rather than starting from ad hoc special cases:
+
+* Match new element constructions against the nearest existing family with the same structural
+  decomposition. Tensor-product, restricted, physically mapped, and enriched elements usually share
+  a factorization pattern that should be reused explicitly.
+* When a feature seems to require a special case, test whether the same mathematics already appears in
+  another element family or mapping path. The right answer is often a more general basis
+  transformation, not a new branch.
+* Separate reference-space reasoning from physical-space reasoning. In FIAT and FInAT, derive basis
+  transformations from the element map and continuity requirements first, then encode that structure
+  in GEM expressions.
+* Treat tensor-product spaces as tensor-product mathematics. Look for Kronecker-style factorization
+  in basis matrices, coordinate mappings, and dual evaluations before introducing custom assembly
+  logic.
+* For extruded or vertically constant factors, identify the dimension that is geometrically active
+  and the dimension that is algebraically passive. The passive factor should usually contribute a
+  simple constant, identity, or lower-dimensional pullback rather than a new geometric rule.
+* Debug by matching the failing object against a known neighboring case: compare the cell, element
+  family, mapping type, continuity class, and tensor structure before changing code.
+* In GEM, inspect whether an expression should factor, broadcast, or propagate a coordinate mapping.
+  If the expression is not matching the expected shape, the bug is often in the way indices or
+  subexpressions are assembled, not in the downstream optimizer.
+* Use the mathematical continuity target as a design constraint. For finite elements, ask what
+  inter-element continuity the space must satisfy, then derive the local basis and transformation
+  rules from that requirement.
+* Prefer proofs by structure over proofs by example. A construction is correct when the pullback,
+  restriction, and tensor-product algebra agree with the element's continuity and approximation
+  properties, not when one or two test cases happen to pass.
+
+## Style and Conventions
+
+When writing Python code for FIAT, maintain the ecosystem's structural and stylistic integrity:
+
+* Class attributes must be declared in one visible location.
+* Attributes must be initialized in the `__init__` constructor or declared as a `functools.cached_property` if they are expensive to compute.
+* Ad hoc lazy initialization discovering attributes via `hasattr`, `setattr`, or `getattr` scattered across methods is strictly prohibited.
+* Boolean attributes must be used to record initialization intent and state instead of probing for the presence of state built by an initialization function.
+* New code must include type hints on all function and method signatures.
+* Public-facing APIs must include properly formatted `numpydoc`-style docstrings.
+* CI enforces `pydocstyle` (see the `[pydocstyle]` section of `setup.cfg` for the active ignore list)
+  in addition to `flake8`; run `pydocstyle <changed files>` locally before finishing a change, since a
+  clean `flake8` pass does not imply a clean `pydocstyle` pass.
diff --git a/FIAT/expansions.py b/FIAT/expansions.py
index 5beebc234..cc517c1ce 100644
--- a/FIAT/expansions.py
+++ b/FIAT/expansions.py
@@ -10,6 +10,7 @@
 import numpy
 import math
 from FIAT import reference_element, jacobi
+from FIAT.precision import calibrate_tolerance
 
 
 def morton_index2(p, q=0):
@@ -62,23 +63,114 @@ def jacobi_factors(x, y, z, dx, dy, dz):
     return fa, fb, fc, dfa, dfb, dfc
 
 
-def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
+def _product_derivative(factor: numpy.ndarray,
+                        dfactor: numpy.ndarray | None,
+                        ddfactor: numpy.ndarray | None,
+                        operands: list[numpy.ndarray],
+                        order: int) -> numpy.ndarray:
+    """Differentiate a recurrence factor times a basis derivative tensor.
+
+    Parameters
+    ----------
+    factor : numpy.ndarray
+        The recurrence factor.
+    dfactor : numpy.ndarray, optional
+        The first derivative of the recurrence factor.
+    ddfactor : numpy.ndarray, optional
+        The second derivative of the recurrence factor.
+    operands : list[numpy.ndarray]
+        List of basis derivative tensors up to order rank.
+    order : int
+        The derivative order of the output.
+
+    Returns
+    -------
+    numpy.ndarray
+        The differentiated product tensor of shape (len(mis(dim, rank)), num_points).
+
+    """
+    from FIAT.polynomial_set import mis
+
+    dim = dfactor.shape[0] if dfactor is not None else 0
+    alphas = mis(dim, order)
+
+    # result = F * D^alpha G
+    result = factor * operands[order]
+    dtype = result.dtype
+
+    if dfactor is not None and order >= 1:
+        alpha_minus1 = mis(dim, order - 1)
+        idx_of_minus1 = {alpha: j for j, alpha in enumerate(alpha_minus1)}
+        DF = numpy.zeros((len(alphas), len(alpha_minus1), *dfactor.shape[1:]), dtype=dtype)
+        for i, alpha in enumerate(alphas):
+            for d in range(dim):
+                if alpha[d] < 1:
+                    continue
+                alpha_minus = list(alpha)
+                alpha_minus[d] -= 1
+                j = idx_of_minus1[tuple(alpha_minus)]
+                DF[i, j] += alpha[d] * dfactor[d]
+        # result += alpha * D F * D^(alpha-1) G
+        result += numpy.einsum("ij...,j...->i...", DF, operands[order-1])
+
+    if ddfactor is not None and order >= 2:
+        alpha_minus2 = mis(dim, order - 2)
+        idx_of_minus2 = {alpha: j for j, alpha in enumerate(alpha_minus2)}
+        DDF = numpy.zeros((len(alphas), len(alpha_minus2), *ddfactor.shape[2:]), dtype=dtype)
+        for i, alpha in enumerate(alphas):
+            for d1 in range(dim):
+                for d2 in range(d1, dim):
+                    if alpha[d1] < 1 + (d1 == d2) or alpha[d2] < 1 + (d1 == d2):
+                        continue
+                    alpha_minus = list(alpha)
+                    alpha_minus[d1] -= 1
+                    alpha_minus[d2] -= 1
+                    j = idx_of_minus2[tuple(alpha_minus)]
+                    if d1 == d2:
+                        a2 = alpha[d1] * (alpha[d1] - 1) // 2
+                    else:
+                        a2 = alpha[d1] * alpha[d2]
+                    DDF[i, j] += a2 * ddfactor[d1, d2]
+        # result += alpha*(alpha-1) * D^2 F * D^(alpha-2) G
+        result += numpy.einsum("ij...,j...->i...", DDF, operands[order-2])
+
+    return result
+
+
+def dubiner_recurrence(dim: int,
+                       n: int,
+                       order: int,
+                       ref_pts: numpy.ndarray,
+                       Jinv: numpy.ndarray,
+                       scale: float,
+                       variant: str | None = None) -> list[numpy.ndarray]:
     """Tabulate a Dubiner expansion set using the recurrence from (Kirby 2010).
 
-    :arg dim: The spatial dimension of the simplex.
-    :arg n: The polynomial degree.
-    :arg order: The maximum order of differentiation.
-    :arg ref_pts: An ``ndarray`` with the coordinates on the default (-1, 1)^d simplex.
-    :arg Jinv: The inverse of the Jacobian of the coordinate mapping from the default simplex.
-    :arg scale: A scale factor that sets the first member of expansion set.
-    :arg variant: Choose between the default (None) orthogonal basis,
-                  'bubble' for integrated Jacobi polynomials,
-                  or 'dual' for the L2-duals of the integrated Jacobi polynomials.
-
-    :returns: A tuple with tabulations of the expansion set and its derivatives.
+    Parameters
+    ----------
+    dim : int
+        The spatial dimension of the simplex.
+    n : int
+        The polynomial degree.
+    order : int
+        The maximum order of differentiation.
+    ref_pts : numpy.ndarray
+        An ``ndarray`` with the coordinates on the default (-1, 1)^d simplex.
+    Jinv : numpy.ndarray
+        The inverse of the Jacobian of the coordinate mapping from the default simplex.
+    scale : float
+        A scale factor that sets the first member of expansion set.
+    variant : str, optional
+        Choose between the default (None) orthogonal basis,
+        'bubble' for integrated Jacobi polynomials,
+        or 'dual' for the L2-duals of the integrated Jacobi polynomials.
+
+    Returns
+    -------
+    list[numpy.ndarray]
+        A list of numpy arrays with tabulations of the expansion set and its derivatives.
+
     """
-    if order > 2:
-        raise ValueError("Higher order derivatives not supported")
     if variant not in [None, "bubble", "dual"]:
         raise ValueError(f"Invalid variant {variant}")
     if variant == "bubble":
@@ -86,16 +178,17 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
 
     num_members = math.comb(n + dim, dim)
 
-    outer = lambda x, y: x[:, None, ...] * y[None, ...]
-
     pad_dim = dim + 2
     dX = pad_jacobian(Jinv, pad_dim)
 
     phi0 = numpy.array([sum((ref_pts[i] - ref_pts[i] for i in range(dim)), 0.0)])
-    results = [numpy.zeros((num_members,) + (len(dX[0]),)*k + phi0.shape[1:], dtype=phi0.dtype)
-               for k in range(order+1)]
+    dtype = phi0.dtype
+    results = [
+        numpy.zeros((num_members, math.comb(len(dX[0])+k-1, k), *phi0.shape[1:]), dtype=dtype)
+        for k in range(order+1)
+    ]
 
-    phi, dphi, ddphi = results + [None] * (2-order)
+    phi = results[0]
     phi[0] = scale
     if dim == 0 or n == 0:
         return results
@@ -109,9 +202,9 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
     for codim in range(dim):
         # Extend the basis from codim to codim + 1
         fa, fb, fc, dfa, dfb, dfc = jacobi_factors(*X[codim:codim+3], *dX[codim:codim+3])
-        ddfc = 2 * outer(dfb, dfb)
+        ddfc = 2 * numpy.outer(dfb, dfb)
         for sub_index in reference_element.lattice_iter(0, n, codim):
-            # handle i = 1
+            # handle i = 0
             icur = idx(*sub_index, 0)
             inext = idx(*sub_index, 1)
 
@@ -126,15 +219,13 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
                 b = 0.5 * (alpha - beta)
 
             fcur = a * fa - b * fb
-            phi[inext] = fcur * phi[icur]
-            if dphi is not None:
+            phi[inext] = phi[icur] * fcur
+            if order:
                 dfcur = a * dfa - b * dfb
-                dphi[inext] = phi[icur] * dfcur
-                dphi[inext] += fcur * dphi[icur]
-                if ddphi is not None:
-                    ddphi[inext] = outer(dphi[icur], dfcur)
-                    ddphi[inext] += outer(dfcur, dphi[icur])
-                    ddphi[inext] += fcur * ddphi[icur]
+                cur = [result[icur] for result in results]
+                deg = sum(sub_index) + 1
+                for k in range(1, min(order, deg)+1):
+                    results[k][inext] = _product_derivative(fcur, dfcur, None, cur, k)
 
             # general i by recurrence
             for i in range(1, n - sum(sub_index)):
@@ -143,28 +234,19 @@ def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale, variant=None):
 
                 fcur = a * fa - b * fb
                 fprev = -c * fc
-                phi[inext] = fcur * phi[icur]
-                phi[inext] += fprev * phi[iprev]
-                if dphi is None:
-                    continue
+                phi[inext] = phi[icur] * fcur
+                phi[inext] += phi[iprev] * fprev
 
                 dfcur = a * dfa - b * dfb
                 dfprev = -c * dfc
-                dphi[inext] = phi[icur] * dfcur
-                dphi[inext] += phi[iprev] * dfprev
-                dphi[inext] += fcur * dphi[icur]
-                dphi[inext] += fprev * dphi[iprev]
-                if ddphi is None:
-                    continue
-
                 ddfprev = -c * ddfc
-                ddphi[inext] = phi[iprev] * ddfprev
-                ddphi[inext] += outer(dphi[icur], dfcur)
-                ddphi[inext] += outer(dfcur, dphi[icur])
-                ddphi[inext] += outer(dphi[iprev], dfprev)
-                ddphi[inext] += outer(dfprev, dphi[iprev])
-                ddphi[inext] += fcur * ddphi[icur]
-                ddphi[inext] += fprev * ddphi[iprev]
+
+                cur = [result[icur] for result in results]
+                prev = [result[iprev] for result in results]
+                deg = sum(sub_index) + 1 + i
+                for k in range(1, min(order, deg)+1):
+                    results[k][inext] = _product_derivative(fcur, dfcur, None, cur, k)
+                    results[k][inext] += _product_derivative(fprev, dfprev, ddfprev, prev, k)
 
         # normalize
         d = codim + 1
@@ -237,7 +319,7 @@ def C0_basis(dim, n, tabulations):
 
         dofs.extend(idx(i, j, k) for k in range(1, n+1) for j in range(1, n-k+1) for i in range(2, n-j-k+1))
 
-    return tuple([phi[i] for i in dofs] for phi in tabulations)
+    return tuple(phi[dofs] for phi in tabulations)
 
 
 def xi_triangle(eta):
@@ -292,7 +374,7 @@ def __init__(self, ref_el, scale=None, variant=None):
         self.scale = scale
         self.variant = variant
         self.continuity = "C0" if variant == "bubble" else None
-        self.recurrence_order = 2
+        self.recurrence_order = math.inf
         self._dmats_cache = {}
         self._cell_node_map_cache = {}
 
@@ -344,18 +426,17 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
             phi = C0_basis(tdim, n, phi)
 
         # Pack linearly independent components into a dictionary
-        result = {(0,) * sd: numpy.asarray(phi[0])}
-        for r in range(1, len(phi)):
-            vr = numpy.transpose(phi[r], tuple(range(1, r+1)) + (0, r+1))
-            for indices in numpy.ndindex(vr.shape[:r]):
-                alpha = tuple(map(indices.count, range(sd)))
-                if alpha not in result:
-                    result[alpha] = vr[indices]
+        result = {}
+        for r in range(len(phi)):
+            vr = numpy.asarray(phi[r])
+            vr = vr.transpose(1, 0, *range(2, vr.ndim))
+            for j, alpha in enumerate(mis(sd, r)):
+                result[alpha] = vr[j]
 
         def distance(alpha, beta):
             return sum(ai != bi for ai, bi in zip(alpha, beta))
 
-        # Only use dmats if tabulate failed
+        # Use dmats only for derivatives above the configured recurrence order.
         for i in range(len(phi), order + 1):
             dmats = self.get_dmats(n, cell=cell)
             for alpha in mis(sd, i):
@@ -380,7 +461,8 @@ def _tabulate(self, n, pts, order=0):
 
         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)
+            tol = calibrate_tolerance(1E-12, numpy.array(self.ref_el.vertices).dtype)
+            Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique, tol=tol)
             for cell, phi in phis.items():
                 for alpha in phi:
                     phi[alpha] *= Xi[cell]
@@ -580,7 +662,7 @@ def __init__(self, ref_el, **kwargs):
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         """Returns a dict of tabulations such that
         tabulations[alpha][i, j] = D^alpha phi_i(pts[j])."""
-        if self.variant is not None:
+        if self.variant is not None or self.ref_el.get_spatial_dimension() != 1:
             return super()._tabulate_on_cell(n, pts, order=order, cell=cell, direction=direction)
 
         A, b = self.affine_mappings[cell]
@@ -738,7 +820,8 @@ def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12):
     :arg ref_el: a SimplicialComplex.
     :arg pt: a physical point on the complex.
     :kwarg unique: Are we assigning a unique cell to points on facets?
-    :kwarg tol: the absolute tolerance.
+    :kwarg tol: the absolute tolerance, already adjusted for the caller's
+        working precision (see `FIAT.precision.calibrate_tolerance`).
     :returns: a list of (weighted) characteristic functions for each subcell.
     """
     import gem
diff --git a/FIAT/hermite.py b/FIAT/hermite.py
index 79c83ff97..15bb4eeee 100644
--- a/FIAT/hermite.py
+++ b/FIAT/hermite.py
@@ -9,60 +9,47 @@
 
 
 class CubicHermiteDualSet(dual_set.DualSet):
-    """The dual basis for Lagrange elements.  This class works for
-    simplices of any dimension.  Nodes are point evaluation at
-    equispaced points."""
-
-    def __init__(self, ref_el):
-        entity_ids = {}
-        nodes = []
-        cur = 0
+    """The dual basis for Hermite elements.  This class works for
+    simplices of any dimension.  Nodes are first order jet at
+    vertices and point evaluation at barycenters of 2D entities."""
 
+    def __init__(self, ref_el, degree, variant=None):
         # make nodes by getting points
         # need to do this dimension-by-dimension, facet-by-facet
         top = ref_el.get_topology()
-        verts = ref_el.get_vertices()
-        sd = ref_el.get_spatial_dimension()
-
-        # get jet at each vertex
+        sd = ref_el.get_topological_dimension()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
-        entity_ids[0] = {}
+        # get first order jet at each vertex
         for v in sorted(top[0]):
-            nodes.append(functional.PointEvaluation(ref_el, verts[v]))
-            pd = functional.PointDerivative
-            for i in range(sd):
-                alpha = [0] * sd
-                alpha[i] = 1
-
-                nodes.append(pd(ref_el, verts[v], alpha))
-
-            entity_ids[0][v] = list(range(cur, cur + 1 + sd))
-            cur += sd + 1
-
-        # now only have dofs at the barycenter, which is the
-        # maximal dimension
-        # no edge dof
-
-        entity_ids[1] = {}
-        for i in top[1]:
-            entity_ids
-            entity_ids[1][i] = []
-
-        if sd > 1:
-            # face dof
-            # point evaluation at barycenter
-            entity_ids[2] = {}
+            pt, = ref_el.make_points(0, v, degree, variant=variant)
+            cur = len(nodes)
+            nodes.append(functional.PointEvaluation(ref_el, pt))
+            if sd == 1:
+                # use normal derivative to support manifolds in 1D
+                nodes.append(functional.PointNormalDerivative(ref_el, v, pt))
+            else:
+                nodes.extend(functional.PointDerivative(ref_el, pt, alpha)
+                             for alpha in polynomial_set.mis(sd, 1))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        if sd == 1:
+            # edge dofs: point evaluations to support higher order in 1D
+            for e in sorted(top[1]):
+                cur = len(nodes)
+                pts = ref_el.make_points(1, e, degree-2, variant=variant)
+                nodes.extend(functional.PointEvaluation(ref_el, pt) for pt in pts)
+                entity_ids[1][e].extend(range(cur, len(nodes)))
+        else:
+            assert degree == 3
+            # no edge dof
+            # face dof: point evaluation at barycenter
             for f in sorted(top[2]):
-                pt = ref_el.make_points(2, f, 3)[0]
-                n = functional.PointEvaluation(ref_el, pt)
-                nodes.append(n)
-                entity_ids[2][f] = list(range(cur, cur + 1))
-                cur += 1
-
-            for dim in range(3, sd + 1):
-                entity_ids[dim] = {}
-                for facet in top[dim]:
-                    entity_ids[dim][facet] = []
+                cur = len(nodes)
+                pt, = ref_el.make_points(2, f, degree, variant=variant)
+                nodes.append(functional.PointEvaluation(ref_el, pt))
+                entity_ids[2][f].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
@@ -70,9 +57,10 @@ def __init__(self, ref_el):
 class CubicHermite(finite_element.CiarletElement):
     """The cubic Hermite finite element.  It is what it is."""
 
-    def __init__(self, ref_el, deg=3):
-        assert deg == 3
-        poly_set = polynomial_set.ONPolynomialSet(ref_el, 3)
-        dual = CubicHermiteDualSet(ref_el)
+    def __init__(self, ref_el, degree=3, variant=None):
+        if variant is None:
+            variant = "gll"
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
+        dual = CubicHermiteDualSet(ref_el, degree, variant=variant)
 
-        super().__init__(poly_set, dual, 3)
+        super().__init__(poly_set, dual, degree)
diff --git a/FIAT/macro.py b/FIAT/macro.py
index 537a1aad1..9ce68f625 100644
--- a/FIAT/macro.py
+++ b/FIAT/macro.py
@@ -487,20 +487,31 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs):
         # Impose C^vorder super-smoothness at interior vertices
         # C^forder automatically gives C^{forder+dim-1} at the interior vertex
         verts = numpy.asarray(ref_el.get_vertices())
+        has_vertex_constraints = False
         for vorder in set(order[0].values()):
             vids = [i for i in order[0] if order[0][i] == vorder]
             facets = chain.from_iterable(ref_el.connectivity[(0, sd-1)][v] for v in vids)
             forder = min(order[sd-1][f] for f in facets)
             sorder = forder + sd - 1
             if vorder > sorder:
+                has_vertex_constraints = True
                 jumps = expansion_set.tabulate_jumps(degree, verts[vids], order=vorder)
                 rows.extend(numpy.vstack(jumps[r].T) for r in range(sorder+1, vorder+1))
 
         if len(rows) > 0:
             for row in rows:
                 row *= 1 / max(numpy.max(abs(row)), 1)
-            dual_mat = numpy.vstack(rows)
-            coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True)
+            if has_vertex_constraints:
+                # Project each constraint block onto the current nullspace so
+                # high-order vertex constraints are not buried in one large SVD.
+                coeffs = numpy.eye(expansion_set.get_num_members(degree))
+                for row in rows:
+                    restricted_row = numpy.dot(row, coeffs.T)
+                    nsp = polynomial_set.spanning_basis(restricted_row, nullspace=True, rtol=1e-12)
+                    coeffs = numpy.dot(nsp, coeffs)
+            else:
+                dual_mat = numpy.vstack(rows)
+                coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True)
         else:
             coeffs = numpy.eye(expansion_set.get_num_members(degree))
 
diff --git a/FIAT/precision.py b/FIAT/precision.py
new file mode 100644
index 000000000..e3bbde626
--- /dev/null
+++ b/FIAT/precision.py
@@ -0,0 +1,20 @@
+"""Tolerance adjustment for the caller's floating-point precision."""
+import math
+
+import numpy
+
+#: Working precision assumed when the caller does not specify one.
+DEFAULT_SCALAR_DTYPE = numpy.dtype("float64")
+
+
+def calibrate_tolerance(tol: float, dtype=DEFAULT_SCALAR_DTYPE) -> float:
+    """Relax `tol` to ``sqrt(tol)`` if `dtype` is single-precision,
+    since float32's unit roundoff (~1.2e-7) is much larger than
+    float64's (~2.2e-16). Otherwise `tol` is returned unchanged.
+
+    :arg tol: the tolerance appropriate for double precision.
+    :arg dtype: the caller's working precision. `None` behaves like
+        double precision.
+    """
+    is_single = numpy.dtype(dtype) == numpy.dtype(numpy.float32)
+    return math.sqrt(tol) if is_single else tol
diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index b896059e1..b91211791 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -94,7 +94,7 @@ def make_lattice(verts, n, interior=0, variant=None):
     family = _decode_family(family)
     D = len(verts)
     X = numpy.array(verts)
-    get_point = lambda alpha: tuple(numpy.dot(_recursive(D - 1, n, alpha, family), X))
+    get_point = lambda alpha: tuple(numpy.dot(_recursive(D - 1, n, alpha, family), X).astype(X.dtype))
     return list(map(get_point, multiindex_equal(D, n, interior)))
 
 
@@ -929,6 +929,17 @@ def get_facet_element(self):
 ReferenceElement = Simplex
 
 
+def cast_vertices(verts, dtype):
+    """Cast a tuple of vertex coordinate tuples to the given numpy dtype.
+
+    :arg verts: a tuple of tuples of vertex coordinates.
+    :arg dtype: a numpy dtype, or `None` to leave `verts` unchanged.
+    """
+    if dtype is None:
+        return verts
+    return tuple(tuple(row) for row in numpy.array(verts, dtype=dtype))
+
+
 class UFCSimplex(Simplex):
 
     def construct_subelement(self, dimension):
@@ -1001,6 +1012,11 @@ def __init__(self):
                     1: edges}
         super().__init__(LINE, verts, topology)
 
+    def compute_normal(self, i):
+        "UFC consistent normal"
+        n = self.compute_tangents(1, 0)[0]
+        return n / numpy.linalg.norm(n)
+
 
 class DefaultTriangle(DefaultSimplex):
     """This is the reference triangle with vertices (-1.0,-1.0),
@@ -1628,19 +1644,27 @@ def make_affine_mapping(xs, ys):
     return AT.T, b
 
 
-def ufc_hypercube(spatial_dim):
+def ufc_hypercube(spatial_dim, dtype=None):
     """Factory function that maps spatial dimension to an instance of
-    the UFC reference hypercube of that dimension."""
+    the UFC reference hypercube of that dimension.
+
+    :arg dtype: optional numpy dtype to cast the vertex coordinates to.
+        Defaults to the plain Python `float` used by the hardcoded
+        vertex coordinates.
+    """
     if spatial_dim == 0:
-        return Point()
+        cell = Point()
     elif spatial_dim == 1:
-        return UFCInterval()
+        cell = UFCInterval()
     elif spatial_dim == 2:
-        return UFCQuadrilateral()
+        cell = UFCQuadrilateral()
     elif spatial_dim == 3:
-        return UFCHexahedron()
+        cell = UFCHexahedron()
     else:
         raise RuntimeError(f"Can't create UFC hypercube of dimension {spatial_dim}.")
+    if dtype is not None:
+        cell.vertices = cast_vertices(cell.vertices, dtype)
+    return cell
 
 
 def default_simplex(spatial_dim):
@@ -1658,34 +1682,42 @@ def default_simplex(spatial_dim):
         raise RuntimeError(f"Can't create default simplex of dimension {spatial_dim}.")
 
 
-def ufc_simplex(spatial_dim):
+def ufc_simplex(spatial_dim, dtype=None):
     """Factory function that maps spatial dimension to an instance of
-    the UFC reference simplex of that dimension."""
+    the UFC reference simplex of that dimension.
+
+    :arg dtype: optional numpy dtype to cast the vertex coordinates to.
+        Defaults to the plain Python `float` used by the hardcoded
+        vertex coordinates.
+    """
     if spatial_dim == 0:
-        return Point()
+        cell = Point()
     elif spatial_dim == 1:
-        return UFCInterval()
+        cell = UFCInterval()
     elif spatial_dim == 2:
-        return UFCTriangle()
+        cell = UFCTriangle()
     elif spatial_dim == 3:
-        return UFCTetrahedron()
+        cell = UFCTetrahedron()
     else:
         raise RuntimeError(f"Can't create UFC simplex of dimension {spatial_dim}.")
+    if dtype is not None:
+        cell.vertices = cast_vertices(cell.vertices, dtype)
+    return cell
 
 
-def symmetric_simplex(spatial_dim):
+def symmetric_simplex(spatial_dim, dtype=None):
     A = numpy.array([[2, 1, 1],
                      [0, numpy.sqrt(3), numpy.sqrt(3)/3],
-                     [0, 0, numpy.sqrt(6)*(2/3)]])
+                     [0, 0, numpy.sqrt(6)*(2/3)]], dtype=dtype)
     A = A[:spatial_dim, :][:, :spatial_dim]
     b = A.sum(axis=1) * (-1 / (1 + spatial_dim))
-    Ref1 = ufc_simplex(spatial_dim)
+    Ref1 = ufc_simplex(spatial_dim, dtype=dtype)
     v = numpy.dot(Ref1.get_vertices(), A.T) + b[None, :]
     vertices = tuple(map(tuple, v))
     return SymmetricSimplex(Ref1.get_shape(), vertices, Ref1.get_topology())
 
 
-def ufc_cell(cell):
+def ufc_cell(cell, dtype=None):
     """Handle incoming calls from FFC."""
 
     # celltype could be a string or a cell.
@@ -1696,19 +1728,19 @@ def ufc_cell(cell):
 
     if " * " in celltype:
         # Tensor product cell
-        return TensorProductCell(*map(ufc_cell, celltype.split(" * ")))
+        return TensorProductCell(*(ufc_cell(c, dtype=dtype) for c in celltype.split(" * ")))
     elif celltype == "quadrilateral":
-        return UFCQuadrilateral()
+        return ufc_hypercube(2, dtype=dtype)
     elif celltype == "hexahedron":
-        return UFCHexahedron()
+        return ufc_hypercube(3, dtype=dtype)
     elif celltype == "vertex":
-        return ufc_simplex(0)
+        return ufc_simplex(0, dtype=dtype)
     elif celltype == "interval":
-        return ufc_simplex(1)
+        return ufc_simplex(1, dtype=dtype)
     elif celltype == "triangle":
-        return ufc_simplex(2)
+        return ufc_simplex(2, dtype=dtype)
     elif celltype == "tetrahedron":
-        return ufc_simplex(3)
+        return ufc_simplex(3, dtype=dtype)
     else:
         raise RuntimeError(f"Don't know how to create UFC cell of type {str(celltype)}")
 
diff --git a/finat/cube.py b/finat/cube.py
index b1517e03f..be3d851bf 100644
--- a/finat/cube.py
+++ b/finat/cube.py
@@ -99,3 +99,6 @@ def value_shape(self):
     @property
     def mapping(self):
         return self.product.mapping
+
+    def dual_evaluation(self, argument, coordinate_mapping=None):
+        return self.product.dual_evaluation(argument, coordinate_mapping)
diff --git a/finat/enriched.py b/finat/enriched.py
index 631dba374..6f5c78f4d 100644
--- a/finat/enriched.py
+++ b/finat/enriched.py
@@ -8,19 +8,29 @@
 from gem.utils import cached_property
 
 from finat.finiteelementbase import FiniteElementBase
+from finat.hdivcurl import HCurlElement, HDivElement
 
 
 class EnrichedElement(FiniteElementBase):
     """A finite element whose basis functions are the union of the
     basis functions of several other finite elements."""
 
-    def __new__(cls, elements):
+    def __new__(cls, elements, is_nodal_enriched=None):
         elements = tuple(chain.from_iterable(e.elements if isinstance(e, EnrichedElement) else (e,) for e in elements))
         if len(elements) == 1:
             return elements[0]
         else:
             self = super().__new__(cls)
             self.elements = elements
+
+            if is_nodal_enriched is None:
+                is_nodal_enriched = all(
+                    is_orthogonal(elements[i], elements[j])
+                    for i in range(len(elements))
+                    for j in range(i+1, len(elements))
+                )
+
+            self.is_nodal_enriched = is_nodal_enriched
             return self
 
     @cached_property
@@ -149,6 +159,34 @@ def mapping(self):
             result, = mappings
             return result
 
+    def dual_evaluation(self, argument, coordinate_mapping=None):
+        if not self.is_nodal_enriched:
+            raise NotImplementedError(
+                f"Dual evaluation not defined for element {type(self).__name__}"
+            )
+        # Gather results from all sub-elements
+        # Each sub_result is (eval_expr, local_indices)
+        sub_results = [sub.dual_evaluation(argument, coordinate_mapping=coordinate_mapping)
+                       for sub in self.elements]
+
+        # Extract the evaluation sub-expressions
+        # We must ensure that all subindices are in the free indices of subexpr
+        # before wrapping in ComponentTensor. If some are missing (e.g. if the
+        # expression simplified to a constant), we multiply by a dummy ones tensor.
+        evals = []
+        for sub, (subexpr, subindices) in zip(self.elements, sub_results):
+            missing_indices = tuple(idx for idx in subindices if idx not in subexpr.free_indices)
+            if missing_indices:
+                shape = tuple(idx.extent for idx in missing_indices)
+                ones = gem.Literal(numpy.ones(shape))
+                dummy = gem.Indexed(ones, missing_indices)
+                subexpr = gem.Product(subexpr, dummy)
+            evals.append(gem.ComponentTensor(subexpr, subindices))
+
+        beta = self.get_indices()
+        expr = gem.Indexed(gem.Concatenate(*evals), beta)
+        return expr, beta
+
 
 def tree_map(f, *args):
     """Like the built-in :py:func:`map`, but applies to a tuple tree."""
@@ -201,3 +239,12 @@ def concatenate_entity_permutations(elements):
                     offset = len(o_e_dim_permutations)
                     o_e_dim_permutations += list(offset + q for q in p)
     return permutations
+
+
+def is_orthogonal(A, B):
+    """Test whether two elements are orthogonal."""
+    if isinstance(A, (HCurlElement, HDivElement)) and isinstance(B, (HCurlElement, HDivElement)):
+        Amap = A.transform(gem.Literal(numpy.ones(A.wrappee.value_shape)))
+        Bmap = B.transform(gem.Literal(numpy.ones(B.wrappee.value_shape)))
+        return sum(a * b for a, b in zip(Amap, Bmap)) == gem.Literal(0.0)
+    return False
diff --git a/finat/hermite.py b/finat/hermite.py
index 5b3a28115..61262c0ab 100644
--- a/finat/hermite.py
+++ b/finat/hermite.py
@@ -1,5 +1,6 @@
 import FIAT
-from gem import ListTensor
+import numpy
+from gem import ListTensor, partial_indexed
 
 from finat.citations import cite
 from finat.fiat_elements import ScalarFiatElement
@@ -7,26 +8,31 @@
 
 
 class Hermite(PhysicallyMappedElement, ScalarFiatElement):
-    def __init__(self, cell, degree=3):
+    def __init__(self, cell, degree=3, variant=None):
         cite("Ciarlet1972")
-        super().__init__(FIAT.CubicHermite(cell))
+        super().__init__(FIAT.CubicHermite(cell, degree=degree, variant=variant))
 
     def basis_transformation(self, coordinate_mapping):
         Js = [coordinate_mapping.jacobian_at(vertex)
               for vertex in self.cell.get_vertices()]
 
+        pns = coordinate_mapping.physical_normals()
         h = coordinate_mapping.cell_size()
 
-        d = self.cell.get_dimension()
         M = identity(self.space_dimension())
 
-        cur = 0
-        for i in range(d+1):
-            cur += 1  # skip the vertex
+        entity_ids = self.entity_dofs()
+        for i in entity_ids[0]:
+            # skip the PointEvaluation DOF
+            vids = entity_ids[0][i][1:]
             J = Js[i]
-            for j in range(d):
-                for k in range(d):
-                    M[cur+j, cur+k] = J[j, k] / h[i]
-            cur += d
+
+            gdim, tdim = J.shape
+            if gdim != tdim:
+                assert tdim == 1
+                J = partial_indexed(pns, (i,)) @ J
+
+            Jnp = numpy.reshape([J[i] for i in numpy.ndindex(J.shape)], J.shape)
+            M[numpy.ix_(vids, vids)] = Jnp * (1 / h[i])
 
         return ListTensor(M)
diff --git a/finat/restricted.py b/finat/restricted.py
index b5144be9f..6b4e291e7 100644
--- a/finat/restricted.py
+++ b/finat/restricted.py
@@ -242,7 +242,7 @@ def restrict_tpe(element, domain, take_closure):
         if all(f is not null_element for f in new_factors):
             elements.append(finat.TensorProductElement(new_factors))
     if elements:
-        return finat.EnrichedElement(elements)
+        return finat.EnrichedElement(elements, is_nodal_enriched=True)
     else:
         return null_element
 
diff --git a/pyproject.toml b/pyproject.toml
index 1a013dac8..e2668d681 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -15,7 +15,7 @@ dependencies = [
   "symengine",
   "sympy",
 ]
-requires-python = ">=3.10"
+requires-python = ">=3.11"
 authors = [
   {name = "Robert C. Kirby et al.", email = "fenics-dev@googlegroups.com"},
   {name = "Imperial College London and others", email = "david.ham@imperial.ac.uk"},
diff --git a/test/FIAT/unit/test_hct.py b/test/FIAT/unit/test_hct.py
index 84e0a2379..cb2d56fb9 100644
--- a/test/FIAT/unit/test_hct.py
+++ b/test/FIAT/unit/test_hct.py
@@ -74,3 +74,9 @@ def test_full_polynomials(cell, reduced):
     C1 = CkPolynomialSet(ref_complex, degree, order=1, variant="bubble")
     C1_tab = C1.tabulate(pts)[(0, 0)]
     assert span_greater_equal(tab, C1_tab)
+
+
+@pytest.mark.parametrize("degree, space_dimension", [(13, 127), (14, 144)])
+def test_hct_high_order_degree(degree: int, space_dimension: int) -> None:
+    fe = HCT(ufc_simplex(2), degree)
+    assert fe.space_dimension() == space_dimension
diff --git a/test/FIAT/unit/test_polynomial.py b/test/FIAT/unit/test_polynomial.py
index c211fb815..efaf50feb 100644
--- a/test/FIAT/unit/test_polynomial.py
+++ b/test/FIAT/unit/test_polynomial.py
@@ -84,6 +84,31 @@ def eval_basis(f, pt):
             assert numpy.allclose(uh, exact, atol=1E-14)
 
 
+@pytest.mark.parametrize("dim", [2, 3])
+@pytest.mark.parametrize("variant", [None, "bubble"])
+def test_high_order_expansion_derivatives(dim, variant):
+    cell = reference_element.default_simplex(dim)
+    degree = 5
+    order = 4
+    points = reference_element.make_lattice(cell.get_vertices(), 5, interior=1)
+
+    fallback = expansions.ExpansionSet(cell, variant=variant)
+    fallback.recurrence_order = 2
+    expected = fallback._tabulate(degree, points, order=order)
+
+    recurrence = expansions.ExpansionSet(cell, variant=variant)
+
+    def get_dmats(*args, **kwargs):
+        raise AssertionError("high-order derivatives should use recurrence tabulation")
+
+    recurrence.get_dmats = get_dmats
+    actual = recurrence._tabulate(degree, points, order=order)
+
+    assert actual.keys() == expected.keys()
+    for alpha in actual:
+        assert numpy.allclose(actual[alpha], expected[alpha], atol=1E-10, rtol=1E-10)
+
+
 @pytest.mark.parametrize("degree", [10])
 def test_expansion_orthonormality(cell, degree):
     U = expansions.ExpansionSet(cell)
diff --git a/test/FIAT/unit/test_precision.py b/test/FIAT/unit/test_precision.py
new file mode 100644
index 000000000..30b82da64
--- /dev/null
+++ b/test/FIAT/unit/test_precision.py
@@ -0,0 +1,72 @@
+import math
+
+import gem
+import numpy
+import pytest
+from gem.interpreter import evaluate
+from gem.node import traversal
+
+from FIAT import DiscontinuousLagrange
+from FIAT.precision import calibrate_tolerance
+from FIAT.reference_element import ufc_simplex
+
+
+@pytest.mark.parametrize("dtype,expected", [
+    (None, 1E-12),
+    (numpy.float64, 1E-12),
+    (numpy.float32, math.sqrt(1E-12)),
+    ("float32", math.sqrt(1E-12)),
+])
+def test_calibrate_tolerance(dtype, expected):
+    assert calibrate_tolerance(1E-12, dtype) == expected
+
+
+def macro_element(dtype):
+    """A macro element whose symbolic tabulation exercises
+    FIAT.expansions.compute_partition_of_unity, where `dtype` is used."""
+    K = ufc_simplex(1, dtype=dtype)
+    return DiscontinuousLagrange(K, 1, variant="iso")
+
+
+def literals(expr):
+    """Return the set of gem.Literal values appearing in a gem expression."""
+    return {node.value for node in traversal([expr]) if isinstance(node, gem.Literal)}
+
+
+def evaluate_at(exprs, x0, value):
+    """Evaluate an iterable of gem expressions at x0 = value."""
+    results = evaluate(list(exprs), bindings={x0: numpy.asarray(value)})
+    return numpy.array([result.arr for result in results])
+
+
+def test_dtype_propagates_into_symbolic_tabulation():
+    """The `dtype`-adjusted tolerance should appear verbatim in the
+    gem expression tree produced by tabulating a macro element,
+    confirming it reaches FIAT.expansions.compute_partition_of_unity."""
+    fe64 = macro_element(dtype=numpy.float64)
+    fe32 = macro_element(dtype=numpy.float32)
+    assert numpy.array(fe64.ref_complex.vertices).dtype == numpy.float64
+    assert numpy.array(fe32.ref_complex.vertices).dtype == numpy.float32
+    x0 = gem.Variable("x0", ())
+
+    tab64 = fe64.tabulate(0, (x0,))[(0,)][0]
+    tab32 = fe32.tabulate(0, (x0,))[(0,)][0]
+
+    assert 1E-12 in literals(tab64)
+    assert math.sqrt(1E-12) in literals(tab32)
+
+    tab64 = fe64.tabulate(0, (x0,))[(0,)]
+    tab32 = fe32.tabulate(0, (x0,))[(0,)]
+
+    # The macro element splits the interval at its midpoint (x0 = 0.5).
+    # 1e-9 lies between the float64 tolerance (1e-12) and the
+    # float32-adjusted tolerance (sqrt(1e-12) = 1e-6), so the two
+    # dtypes classify this point into different subcells.
+    near_boundary = 0.5 + 1e-9
+    assert not numpy.allclose(evaluate_at(tab64, x0, near_boundary),
+                              evaluate_at(tab32, x0, near_boundary))
+
+    # Far from the boundary, both dtypes classify the point the same way.
+    away_from_boundary = 0.5 + 1e-3
+    assert numpy.allclose(evaluate_at(tab64, x0, away_from_boundary),
+                          evaluate_at(tab32, x0, away_from_boundary))
diff --git a/test/finat/conftest.py b/test/finat/conftest.py
index 1fab8dd3d..a4d875253 100644
--- a/test/finat/conftest.py
+++ b/test/finat/conftest.py
@@ -107,18 +107,19 @@ def scaled_simplex(dim, scale):
 
 @pytest.fixture
 def ref_el():
-    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
+    K = {dim: FIAT.ufc_simplex(dim) for dim in (1, 2, 3)}
     return K
 
 
 @pytest.fixture
 def phys_el():
-    K = {dim: FIAT.ufc_simplex(int(dim)) for dim in (2, 2.5, 3)}
+    K = {dim: FIAT.ufc_simplex(int(dim)) for dim in (1.5, 2, 2.5, 3)}
     K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))
     K[3].vertices = ((0, 0, 0),
                      (1., 0.1, -0.37),
                      (0.01, 0.987, -.23),
                      (-0.1, -0.2, 1.38))
+    K[1.5].vertices = K[2].vertices[:-1]
     K[2.5].vertices = K[3].vertices[:-1]
     return K
 
diff --git a/test/finat/test_dual_basis.py b/test/finat/test_dual_basis.py
index dba58251c..7cfe1ae6a 100644
--- a/test/finat/test_dual_basis.py
+++ b/test/finat/test_dual_basis.py
@@ -1,6 +1,7 @@
 import pytest
 import numpy
 import finat
+import gem
 from FIAT import ufc_simplex
 
 
@@ -27,3 +28,20 @@ def test_collapse_repeated_points(dim):
 
     assert len(points) == len(numpy.unique(numpy.round(points, decimals=7), axis=0))
     assert len(points) == expected
+
+
+def test_enriched_element_dual_evaluation():
+    cell = ufc_simplex(2)
+    fe = finat.Lagrange(cell, 3)
+
+    fe1 = finat.RestrictedElement(fe, restriction_domain="interior")
+    fe2 = finat.RestrictedElement(fe, restriction_domain="facet")
+    enriched = finat.EnrichedElement([fe1, fe2], is_nodal_enriched=True)
+
+    # Check that calling dual_evaluation returns a valid Indexed expression
+    fn = lambda x: gem.Literal(1.0)
+    expr, indices = enriched.dual_evaluation(fn)
+    assert isinstance(expr, gem.Indexed)
+    assert isinstance(expr.children[0], gem.Concatenate)
+    assert len(indices) == 1
+    assert indices[0].extent == enriched.space_dimension()
diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py
index fe74d2ad4..eb22eb82a 100644
--- a/test/finat/test_zany_mapping.py
+++ b/test/finat/test_zany_mapping.py
@@ -119,6 +119,13 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs):
     assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), pp.pformat((np.round(error, 8).tolist(), *inds))
 
 
+@pytest.mark.parametrize("element, degree", [
+    *((finat.Hermite, k) for k in range(3, 6)),
+])
+def test_C1_interval(ref_to_phys, element, degree):
+    check_zany_mapping(element, ref_to_phys[1.5], degree)
+
+
 @pytest.mark.parametrize("element", [
                          finat.Morley,
                          finat.Hermite,
@@ -132,6 +139,7 @@ def test_C1_triangle(ref_to_phys, element):
 
 @pytest.mark.parametrize("element", [
                          finat.Morley,
+                         finat.Hermite,
                          finat.Walkington,
                          ])
 def test_C1_tetrahedron(ref_to_phys, element):