firedrakeproject · pbrubeck · Jan 21, 2026 · Jan 13, 2026 · Jan 13, 2026 · Jan 13, 2026
diff --git a/FIAT/__init__.py b/FIAT/__init__.py
@@ -62,6 +62,7 @@
 from FIAT.bubble import Bubble, FacetBubble
 from FIAT.hdiv_trace import HDivTrace
 from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen
+from FIAT.wuxu import WuXuH3NC, WuXuRobustH3NC
 from FIAT.walkington import Walkington
 from FIAT.histopolation import Histopolation
 from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite  # noqa: F401
@@ -88,6 +89,8 @@
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
                       "Hermite": CubicHermite,
+                      "Nonconforming Wu-Xu": WuXuH3NC,
+                      "Nonconforming Robust Wu-Xu": WuXuRobustH3NC,
                       "Hsieh-Clough-Tocher": HsiehCloughTocher,
                       "QuadraticPowellSabin6": QuadraticPowellSabin6,
                       "QuadraticPowellSabin12": QuadraticPowellSabin12,

diff --git a/FIAT/argyris.py b/FIAT/argyris.py
@@ -7,8 +7,7 @@
 from FIAT import finite_element, polynomial_set, dual_set
 from FIAT.check_format_variant import check_format_variant, parse_quadrature_scheme
 from FIAT.functional import (PointEvaluation, PointDerivative, PointNormalDerivative,
-                             IntegralMoment,
-                             IntegralMomentOfNormalDerivative)
+                             IntegralMoment, IntegralMomentOfDerivative)
 from FIAT.jacobi import eval_jacobi_batch, eval_jacobi_deriv_batch
 from FIAT.quadrature import FacetQuadratureRule
 from FIAT.reference_element import TRIANGLE, ufc_simplex
@@ -37,16 +36,17 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
             # edge dofs
             k = degree - 5
             rline = ufc_simplex(1)
-            Q = parse_quadrature_scheme(rline, interpolant_deg+k-1, quad_scheme)
-            x = 2.0 * Q.get_points() - 1.0
-            phis = eval_jacobi_batch(2, 2, k, x)
-            dphis = eval_jacobi_deriv_batch(2, 2, k, x)
+            Q_ref = parse_quadrature_scheme(rline, interpolant_deg+k-1, quad_scheme)
+            x = rline.compute_barycentric_coordinates(Q_ref.get_points())
+            xref = x[:, [1]] - x[:, [0]]
+            phis = eval_jacobi_batch(2, 2, k, xref)
+            dphis = 2*eval_jacobi_deriv_batch(2, 2, k, xref)
             for e in sorted(top[1]):
-                Q_mapped = FacetQuadratureRule(ref_el, 1, e, Q)
-                scale = 2 / Q_mapped.jacobian_determinant()
+                Q = FacetQuadratureRule(ref_el, 1, e, Q_ref, avg=True)
+                n = ref_el.compute_normal(e)
                 cur = len(nodes)
-                nodes.extend(IntegralMomentOfNormalDerivative(ref_el, e, Q, phi) for phi in phis)
-                nodes.extend(IntegralMoment(ref_el, Q_mapped, dphi * scale) for dphi in dphis[1:])
+                nodes.extend(IntegralMomentOfDerivative(ref_el, Q, phi, n) for phi in phis)
+                nodes.extend(IntegralMoment(ref_el, Q, dphi) for dphi in dphis[1:])
                 entity_ids[1][e].extend(range(cur, len(nodes)))
 
             # interior dofs
@@ -55,11 +55,9 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
                 cell = ref_el.construct_subelement(sd)
                 Q_ref = parse_quadrature_scheme(cell, interpolant_deg + q, quad_scheme)
                 Pq = polynomial_set.ONPolynomialSet(cell, q, scale=1)
-                Phis = Pq.tabulate(Q_ref.get_points())[(0,) * sd]
+                phis = Pq.tabulate(Q_ref.get_points())[(0,) * sd]
                 for entity in sorted(top[sd]):
-                    Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
-                    scale = 1 / Q.jacobian_determinant()
-                    phis = scale * Phis
+                    Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref, avg=True)
                     cur = len(nodes)
                     nodes.extend(IntegralMoment(ref_el, Q, phi) for phi in phis)
                     entity_ids[sd][entity] = list(range(cur, len(nodes)))

diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py
@@ -34,9 +34,8 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
             Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)]
             for f in top[sd - 1]:
                 cur = len(nodes)
-                Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref)
-                Jdet = Q.jacobian_determinant()
-                n = ref_el.compute_scaled_normal(f) / Jdet
+                Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref, avg=True)
+                n = ref_el.compute_scaled_normal(f)
                 phis = n[None, :, None] * Pq_at_qpts[:, None, :]
                 nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi)
                              for phi in phis)

diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py
@@ -38,20 +38,18 @@ def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
                 facet = ref_el.construct_subelement(dim)
                 if dim == 0:
                     Q_facet = parse_quadrature_scheme(facet, degree + interpolant_deg-1, quad_scheme)
-                    Phis = numpy.ones((1, len(Q_facet.pts)))
+                    phis = numpy.ones((1, len(Q_facet.pts)))
                 else:
                     k = degree - 1 if dim == sd-1 else degree - (1+dim)
                     if k < 0:
                         continue
                     Q_facet = parse_quadrature_scheme(facet, k + interpolant_deg, quad_scheme)
                     poly_set = polynomial_set.ONPolynomialSet(facet, k)
-                    Phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim]
+                    phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim]
 
                 for i in sorted(top[dim]):
                     cur = len(nodes)
-                    Q = FacetQuadratureRule(ref_el, dim, i, Q_facet)
-                    scale = 1 / Q.jacobian_determinant()
-                    phis = scale * Phis
+                    Q = FacetQuadratureRule(ref_el, dim, i, Q_facet, avg=True)
                     nodes.extend(functional.IntegralMoment(ref_el, Q, phi) for phi in phis)
                     entity_ids[dim][i].extend(range(cur, len(nodes)))
         else:

diff --git a/FIAT/functional.py b/FIAT/functional.py
@@ -13,9 +13,10 @@
 # - type information
 
 from itertools import chain
+from collections import defaultdict
 import numpy
 
-from FIAT import polynomial_set, jacobi, quadrature_schemes
+from FIAT import polynomial_set, jacobi, quadrature, quadrature_schemes
 
 
 class Functional(object):
@@ -242,21 +243,16 @@ def __init__(self, ref_el, edge_no, pt, comp=(), shp=()):
 class PointSecondDerivative(Functional):
     """Represents d/ds1 d/ds2 at a point."""
     def __init__(self, ref_el, s1, s2, pt, comp=(), shp=(), nm=None):
+        S = numpy.outer(s1, s2)
         sd = ref_el.get_spatial_dimension()
-        tau = numpy.zeros((sd*(sd+1)//2,))
-
-        alphas = []
-        cur = 0
-        for i in range(sd):
-            for j in range(i, sd):
-                alpha = [0] * sd
+        tau = defaultdict(float)
+        for index in numpy.ndindex(S.shape):
+            alpha = [0] * sd
+            for i in index:
                 alpha[i] += 1
-                alpha[j] += 1
-                alphas.append(tuple(alpha))
-                tau[cur] = s1[i] * s2[j] + (i != j) * s2[i] * s1[j]
-                cur += 1
+            tau[tuple(alpha)] += S[index]
 
-        dpt_dict = {tuple(pt): [(tau[i], alphas[i], comp) for i in range(len(alphas))]}
+        dpt_dict = {tuple(pt): [(tau[alpha], alpha, comp) for alpha in tau]}
 
         super().__init__(ref_el, shp, {}, dpt_dict, nm or "PointSecondDeriv")
 
@@ -320,82 +316,52 @@ def __call__(self, fn):
 
 
 class IntegralMomentOfDerivative(Functional):
-    """Functional giving directional derivative integrated against some function on a facet."""
+    """Functional giving multi-directional derivative integrated against some function.
 
-    def __init__(self, ref_el, s, Q, f_at_qpts, comp=(), shp=()):
-        self.f_at_qpts = f_at_qpts
+    :arg ref_el: a :class:`Cell`.
+    :arg Q: a :class:`QuadratureRule`.
+    :arg f_at_qpts: an array tabulating the function f at the quadrature
+         points.
+    :arg *directions: a list of vectors of directions of differentiation.
+    :arg comp: Optional argument indicating that only a particular
+         component of the input function should be integrated against f
+    :arg shp: Optional argument giving the value shape of input functions.
+    """
+
+    def __init__(self, ref_el, Q, f_at_qpts, *directions, comp=(), shp=(), nm=""):
         self.Q = Q
+        self.f_at_qpts = f_at_qpts
+        self.comp = comp
+
+        S = directions[0]
+        for dj in directions[1:]:
+            S = numpy.outer(S, dj)
 
         sd = ref_el.get_spatial_dimension()
+        tau = defaultdict(float)
+        for index in numpy.ndindex(S.shape):
+            alpha = [0] * sd
+            for i in index:
+                alpha[i] += 1
+            tau[tuple(alpha)] += S[index]
 
         points = Q.get_points()
         weights = numpy.multiply(f_at_qpts, Q.get_weights())
-
-        alphas = tuple(map(tuple, numpy.eye(sd, dtype=int)))
-        dpt_dict = {tuple(pt): [(wt*s[i], alphas[i], comp) for i in range(sd)]
+        dpt_dict = {tuple(pt): [(wt*tau[alpha], alpha, comp) for alpha in tau]
                     for pt, wt in zip(points, weights)}
 
-        super().__init__(ref_el, shp,
-                         {}, dpt_dict, "IntegralMomentOfDerivative")
+        super().__init__(ref_el, shp, {}, dpt_dict, nm or "IntegralMomentOfDerivative")
 
 
-class IntegralMomentOfNormalDerivative(Functional):
+class IntegralMomentOfNormalDerivative(IntegralMomentOfDerivative):
     """Functional giving normal derivative integrated against some function on a facet."""
 
-    def __init__(self, ref_el, facet_no, Q, f_at_qpts, n=None):
-        if n is None:
-            n = ref_el.compute_normal(facet_no)
-        self.n = n
-        self.f_at_qpts = f_at_qpts
-        self.Q = Q
-
-        sd = ref_el.get_spatial_dimension()
-
+    def __init__(self, ref_el, facet_no, Q_face, f_at_qpts):
+        n = ref_el.compute_normal(facet_no)
         # map points onto facet
-        transform = ref_el.get_entity_transform(sd-1, facet_no)
-        points = transform(Q.get_points())
-        self.dpts = points
-        weights = numpy.multiply(f_at_qpts, Q.get_weights())
-
-        alphas = tuple(map(tuple, numpy.eye(sd, dtype=int)))
-        dpt_dict = {tuple(pt): [(wt*n[i], alphas[i], tuple()) for i in range(sd)]
-                    for pt, wt in zip(points, weights)}
-
-        super().__init__(ref_el, tuple(),
-                         {}, dpt_dict, "IntegralMomentOfNormalDerivative")
-
-
-class IntegralMomentOfBidirectionalDerivative(Functional):
-    """Functional giving second derivative integrated against some function on a facet."""
-
-    def __init__(self, ref_el, Q, f_at_qpts, s1, s2):
-        self.f_at_qpts = f_at_qpts
-        self.Q = Q
-
         sd = ref_el.get_spatial_dimension()
-
-        # map points onto facet
-        points = Q.get_points()
-        self.dpts = points
-        weights = numpy.multiply(f_at_qpts, Q.get_weights())
-
-        tau = numpy.zeros((sd*(sd+1)//2,))
-        alphas = []
-        cur = 0
-        for i in range(sd):
-            for j in range(i, sd):
-                alpha = [0] * sd
-                alpha[i] += 1
-                alpha[j] += 1
-                alphas.append(tuple(alpha))
-                tau[cur] = s1[i] * s2[j] + (i != j) * s2[i] * s1[j]
-                cur += 1
-
-        dpt_dict = {tuple(pt): [(wt*tau[i], alphas[i], tuple()) for i in range(len(alphas))]
-                    for pt, wt in zip(points, weights)}
-
-        super().__init__(ref_el, tuple(),
-                         {}, dpt_dict, "IntegralMomentOfBidirectionalDerivative")
+        Q = quadrature.FacetQuadratureRule(ref_el, sd-1, facet_no, Q_face, avg=True)
+        super().__init__(ref_el, Q, f_at_qpts, n, nm="IntegralMomentOfNormalDerivative")
 
 
 class FrobeniusIntegralMoment(IntegralMoment):

diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py
@@ -26,13 +26,11 @@ def __init__(self, ref_el, degree, quad_scheme=None):
 
             for entity in sorted(top[dim]):
                 cur = len(nodes)
-                Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-                Jdet = Q.jacobian_determinant()
+                Q = FacetQuadratureRule(ref_el, dim, entity, Q_ref, avg=True)
                 for f in ref_el.get_connectivity()[(dim, sd-1)][entity]:
                     normal = ref_el.compute_scaled_normal(f)
                     tangents = ref_el.compute_tangents(sd-1, f)
-                    n = normal / Jdet
-                    nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi)
+                    nodes.extend(BidirectionalMoment(ref_el, t, normal, Q, phi)
                                  for phi in phis for t in tangents)
                 entity_ids[dim][entity].extend(range(cur, len(nodes)))
 

diff --git a/FIAT/hct.py b/FIAT/hct.py
@@ -7,7 +7,7 @@
 # Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024
 
 from FIAT.functional import (PointEvaluation, PointDerivative,
-                             IntegralMoment,
+                             IntegralMoment, IntegralMomentOfDerivative,
                              IntegralMomentOfNormalDerivative)
 from FIAT import finite_element, dual_set, macro, polynomial_set
 from FIAT.check_format_variant import parse_quadrature_scheme
@@ -42,34 +42,34 @@ def __init__(self, ref_complex, degree, reduced=False, quad_scheme=None):
 
         k = 2 if reduced else degree - 3
         facet = ufc_simplex(1)
-        Q = parse_quadrature_scheme(facet, degree-1+k, quad_scheme)
-        qpts = Q.get_points()
-        xref = 2.0 * qpts - 1.0
+        Q_ref = parse_quadrature_scheme(facet, degree-1+k, quad_scheme)
+        x = facet.compute_barycentric_coordinates(Q_ref.get_points())
+        xref = x[:, [1]] - x[:, [0]]
         if reduced:
             f_at_qpts = eval_jacobi(0, 0, k, xref[:, 0])
             for e in sorted(top[1]):
                 cur = len(nodes)
-                nodes.append(IntegralMomentOfNormalDerivative(ref_el, e, Q, f_at_qpts))
+                nodes.append(IntegralMomentOfNormalDerivative(ref_el, e, Q_ref, f_at_qpts))
                 entity_ids[1][e].extend(range(cur, len(nodes)))
         else:
             phis = eval_jacobi_batch(1, 1, k, xref)
-            dphis = eval_jacobi_deriv_batch(1, 1, k, xref)
+            dphis = 2*eval_jacobi_deriv_batch(1, 1, k, xref)
             for e in sorted(top[1]):
-                Q_mapped = FacetQuadratureRule(ref_el, 1, e, Q)
-                scale = 2 / Q_mapped.jacobian_determinant()
+                Q = FacetQuadratureRule(ref_el, 1, e, Q_ref, avg=True)
+                n = ref_el.compute_normal(e)
                 cur = len(nodes)
-                nodes.extend(IntegralMomentOfNormalDerivative(ref_el, e, Q, phi) for phi in phis)
-                nodes.extend(IntegralMoment(ref_el, Q_mapped, dphi * scale) for dphi in dphis[1:])
+                nodes.extend(IntegralMomentOfDerivative(ref_el, Q, phi, n) for phi in phis)
+                nodes.extend(IntegralMoment(ref_el, Q, dphi) for dphi in dphis[1:])
                 entity_ids[1][e].extend(range(cur, len(nodes)))
 
             q = degree - 4
             if q >= 0:
                 Q = parse_quadrature_scheme(ref_complex, degree + q, quad_scheme)
                 Pq = polynomial_set.ONPolynomialSet(ref_el, q, scale=1)
                 phis = Pq.tabulate(Q.get_points())[(0,) * sd]
-                scale = 1 / ref_el.volume()
+                phis *= 1/ref_el.volume()
                 cur = len(nodes)
-                nodes.extend(IntegralMoment(ref_el, Q, phi * scale) for phi in phis)
+                nodes.extend(IntegralMoment(ref_el, Q, phi) for phi in phis)
                 entity_ids[sd][0] = list(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)

diff --git a/FIAT/hellan_herrmann_johnson.py b/FIAT/hellan_herrmann_johnson.py
@@ -67,10 +67,9 @@ def __init__(self, ref_el, degree, variant, qdegree, quad_scheme):
 
             for f in sorted(top[sd-1]):
                 cur = len(nodes)
-                Q = FacetQuadratureRule(ref_el, sd-1, f, Q_ref)
-                detJ = Q.jacobian_determinant()
+                Q = FacetQuadratureRule(ref_el, sd-1, f, Q_ref, avg=True)
                 # n[f]^T u n[f] integrated against a basis for Pk
-                nodes.extend(BidirectionalMoment(ref_el, n[f], n[f]/detJ, Q, phi) for phi in Phis)
+                nodes.extend(BidirectionalMoment(ref_el, n[f], n[f], Q, phi) for phi in Phis)
                 entity_ids[sd-1][f].extend(range(cur, len(nodes)))
 
             ref_facet = ref_el.construct_subelement(sd)
@@ -83,13 +82,12 @@ def __init__(self, ref_el, degree, variant, qdegree, quad_scheme):
             for entity in sorted(top[sd]):
                 cur = len(nodes)
                 faces = cell_to_faces[entity]
-                Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref)
-                detJ = Q.jacobian_determinant()
+                Q = FacetQuadratureRule(ref_el, sd, entity, Q_ref, avg=True)
                 # n[f]^T u n[f] integrated against a basis for P_{k-1}
-                nodes.extend(BidirectionalMoment(ref_el, n[f], n[f]/detJ, Q, phi)
+                nodes.extend(BidirectionalMoment(ref_el, n[f], n[f], Q, phi)
                              for phi in Phis[:dimPkm1] for f in faces)
                 # n[i+1]^T u n[i+2] integrated against a basis for Pk
-                nodes.extend(BidirectionalMoment(ref_el, n[faces[i+1]], n[faces[i+2]]/detJ, Q, phi)
+                nodes.extend(BidirectionalMoment(ref_el, n[faces[i+1]], n[faces[i+2]], Q, phi)
                              for phi in Phis for i in range((sd-1)*(sd-2)))
                 entity_ids[sd][entity].extend(range(cur, len(nodes)))
 

diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py
@@ -47,13 +47,10 @@ def __init__(self, ref_el, degree, codim=0, interpolant_deg=None, quad_scheme=No
         ref_facet = ref_el.construct_subelement(dim)
         poly_set = ONPolynomialSet(ref_facet, degree, scale="L2 piola")
         Q_ref = parse_quadrature_scheme(ref_facet, degree + interpolant_deg, quad_scheme)
-        Phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim]
+        phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim]
         for entity in sorted(top[dim]):
             cur = len(nodes)
-            Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-            # phis must transform like a d-form to undo the measure transformation
-            scale = 1 / Q_facet.jacobian_determinant()
-            phis = scale * Phis
+            Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref, avg=True)
             nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
             entity_ids[dim][entity].extend(range(cur, len(nodes)))
 
@@ -93,13 +90,10 @@ def __init__(self, ref_el, degree, interpolant_deg=None, quad_scheme=None):
             if degree <= dim:
                 continue
             ref_facet = symmetric_simplex(dim)
-            Q_ref, Phis = make_dual_bubbles(ref_facet, degree, interpolant_deg=interpolant_deg, quad_scheme=quad_scheme)
+            Q_ref, phis = make_dual_bubbles(ref_facet, degree, interpolant_deg=interpolant_deg, quad_scheme=quad_scheme)
             for entity in sorted(top[dim]):
                 cur = len(nodes)
-                Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref)
-                # phis must transform like a d-form to undo the measure transformation
-                scale = 1 / Q_facet.jacobian_determinant()
-                phis = scale * Phis
+                Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref, avg=True)
                 nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis)
                 entity_ids[dim][entity].extend(range(cur, len(nodes)))