C2 (macro)elements

FIAT/__init__.py

@@ -22,6 +22,7 @@
 from FIAT.bernstein import Bernstein
 from FIAT.bell import Bell
 from FIAT.hct import HsiehCloughTocher
+from FIAT.double_alfeld import DoubleAlfeld
 from FIAT.alfeld_sorokina import AlfeldSorokina
 from FIAT.arnold_qin import ArnoldQin
 from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div
@@ -62,6 +63,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,9 +90,12 @@
                       "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,
+                      "Double Alfeld": DoubleAlfeld,
                       "Alfeld-Sorokina": AlfeldSorokina,
                       "Arnold-Qin": ArnoldQin,
                       "Christiansen-Hu": ChristiansenHu,

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)))

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)

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:

FIAT/double_alfeld.py

@@ -0,0 +1,94 @@
+# Copyright (C) 2026 Pablo D. Brubeck
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+#
+# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2026
+
+from FIAT.functional import (PointEvaluation, PointDerivative,
+                             IntegralMoment, IntegralMomentOfDerivative)
+from FIAT import finite_element, dual_set, macro, polynomial_set
+from FIAT.polynomial_set import mis
+from FIAT.check_format_variant import parse_quadrature_scheme
+from FIAT.reference_element import TRIANGLE, ufc_simplex
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.jacobi import eval_jacobi_batch, eval_jacobi_deriv_batch
+
+
+class DoubleAlfeldDualSet(dual_set.DualSet):
+    def __init__(self, ref_complex, degree, reduced=False, quad_scheme=None):
+        if degree < 5:
+            raise ValueError("Double Alfeld only defined for degree >= 5")
+        ref_el = ref_complex.get_parent()
+        if ref_el.get_shape() != TRIANGLE:
+            raise ValueError("Double Alfeld only defined on triangles")
+        top = ref_el.get_topology()
+        verts = ref_el.get_vertices()
+        sd = ref_el.get_spatial_dimension()
+        entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
+
+        # get second order jet at each vertex
+        order = 2
+        nodes = []
+        for v in sorted(top[0]):
+            pt = verts[v]
+            cur = len(nodes)
+            nodes.append(PointEvaluation(ref_el, pt))
+            nodes.extend(PointDerivative(ref_el, pt, alpha) for i in range(1, order+1) for alpha in mis(sd, i))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        k = degree-1 if reduced else degree-4
+        facet = ufc_simplex(1)
+        Q_ref = parse_quadrature_scheme(facet, degree-2+k, quad_scheme)
+        x = facet.compute_barycentric_coordinates(Q_ref.get_points())
+        xref = x[:, [1]] - x[:, [0]]
+
+        if reduced:
+            raise NotImplementedError
+        else:
+            phis = eval_jacobi_batch(order, order, k, xref)
+            dphis = 2*eval_jacobi_deriv_batch(order, order, k, xref, order=1)
+            ddphis = 4*eval_jacobi_deriv_batch(order, order, k, xref, order=2)
+            for e in sorted(top[1]):
+                Q = FacetQuadratureRule(ref_el, 1, e, Q_ref, avg=True)
+                n = ref_el.compute_normal(e)
+                cur = len(nodes)
+                nodes.extend(IntegralMoment(ref_el, Q, ddphi) for ddphi in ddphis[2:])
+                nodes.extend(IntegralMomentOfDerivative(ref_el, Q, dphi, n) for dphi in dphis[1:])
+                nodes.extend(IntegralMomentOfDerivative(ref_el, Q, phi, n, n) for phi in phis)
+                entity_ids[1][e].extend(range(cur, len(nodes)))
+
+            q = degree - 6
+            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]
+                phis *= 1/ref_el.volume()
+                cur = len(nodes)
+                nodes.extend(IntegralMoment(ref_el, Q, phi) for phi in phis)
+                entity_ids[sd][0].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class DoubleAlfeld(finite_element.CiarletElement):
+    """The double Alfeld C^2 macroelement on a double barycentric split.
+    See Section 7.5 of Lai & Schumacher for the quintic C2 spline.
+    """
+    def __init__(self, ref_el, degree=5, reduced=False, quad_scheme=None):
+        # Construct the quintic C2 spline on the double Alfeld split
+        ref_complex = macro.AlfeldSplit(macro.AlfeldSplit(ref_el))
+        dual = DoubleAlfeldDualSet(ref_complex, degree, reduced=reduced, quad_scheme=quad_scheme)
+
+        # C3 on major split facets, C2 elsewhere
+        order = {}
+        order[1] = dict.fromkeys(ref_complex.get_interior_facets(1), 2)
+        order[1].update(dict.fromkeys(range(3, 6), degree-2))
+        # C4 at minor split barycenters, C3 at major split barycenter
+        order[0] = dict.fromkeys(ref_complex.get_interior_facets(0), degree-1)
+        order[0][3] = degree-2
+
+        poly_set = macro.CkPolynomialSet(ref_complex, degree, order=order, variant="bubble")
+        formdegree = 0
+        super().__init__(poly_set, dual, degree, formdegree=formdegree)

FIAT/dual_set.py

@@ -7,6 +7,8 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 import numpy
+from itertools import chain
+from collections import defaultdict
 
 from FIAT import polynomial_set, functional
 from FIAT.reference_element import compute_unflattening_map
@@ -117,36 +119,37 @@ def to_riesz(self, poly_set):
         riesz_shape = (num_nodes, *tshape, num_exp)
         mat = numpy.zeros(riesz_shape, "d")
 
-        pts = set()
-        dpts = set()
-        Qs_to_ells = dict()
-        for i, ell in enumerate(self.nodes):
-            if len(ell.deriv_dict) > 0:
-                dpts.update(ell.deriv_dict.keys())
-                continue
-            if isinstance(ell, functional.IntegralMoment):
-                Q = ell.Q
-            else:
-                Q = None
-                pts.update(ell.pt_dict.keys())
-            if Q in Qs_to_ells:
+        def map_quadratures_to_points(nodes, deriv=False):
+            Qs_to_ells = defaultdict(list)
+            for i, ell in enumerate(nodes):
+                if deriv and len(ell.deriv_dict) == 0:
+                    continue
+                elif not deriv and len(ell.pt_dict) == 0:
+                    continue
+                if isinstance(ell, (functional.IntegralMoment, functional.IntegralMomentOfDerivative)):
+                    Q = ell.Q
+                else:
+                    Q = None
                 Qs_to_ells[Q].append(i)
-            else:
-                Qs_to_ells[Q] = [i]
-
-        Qs_to_pts = {}
-        if len(pts) > 0:
-            Qs_to_pts[None] = tuple(sorted(pts))
-        for Q in Qs_to_ells:
-            if Q is not None:
-                cur_pts = tuple(map(tuple, Q.pts))
+            pts = set()
+            Qs_to_pts = {}
+            for Q in Qs_to_ells:
+                if Q is None:
+                    if deriv:
+                        cur_pts = chain.from_iterable(nodes[i].deriv_dict.keys() for i in Qs_to_ells[None])
+                    else:
+                        cur_pts = chain.from_iterable(nodes[i].pt_dict.keys() for i in Qs_to_ells[None])
+                    cur_pts = tuple(set(cur_pts))
+                else:
+                    cur_pts = tuple(map(tuple, Q.pts))
                 Qs_to_pts[Q] = cur_pts
                 pts.update(cur_pts)
+            pts = list(sorted(pts))
+            return Qs_to_ells, Qs_to_pts, pts
 
         # Now tabulate the function values
-        pts = list(sorted(pts))
+        Qs_to_ells, Qs_to_pts, pts = map_quadratures_to_points(self.nodes)
         expansion_values = numpy.transpose(es.tabulate(ed, pts))
-
         for Q in Qs_to_ells:
             ells = Qs_to_ells[Q]
             cur_pts = Qs_to_pts[Q]
@@ -171,25 +174,35 @@ def to_riesz(self, poly_set):
         # Tabulate the derivative values that are needed
         max_deriv_order = max(ell.max_deriv_order for ell in self.nodes)
         if max_deriv_order > 0:
-            dpts = list(sorted(dpts))
+            Qs_to_ells, Qs_to_pts, pts = map_quadratures_to_points(self.nodes, deriv=True)
             # It's easiest/most efficient to get derivatives of the
             # expansion set through the polynomial set interface.
             # This is creating a short-lived set to do just this.
             coeffs = numpy.eye(num_exp)
             expansion = polynomial_set.PolynomialSet(self.ref_el, ed, ed, es, coeffs)
-            dexpansion_values = expansion.tabulate(dpts, max_deriv_order)
-
-            ells = [k for k, ell in enumerate(self.nodes) if len(ell.deriv_dict) > 0]
-            wshape = (len(ells), *tshape, len(dpts))
-            dwts = {alpha: numpy.zeros(wshape, "d") for alpha in dexpansion_values if sum(alpha) > 0}
-            for i, k in enumerate(ells):
-                ell = self.nodes[k]
-                for pt, wac_list in ell.deriv_dict.items():
-                    j = dpts.index(pt)
-                    for (w, alpha, c) in wac_list:
-                        dwts[alpha][i][c][j] = w
-            for alpha in dwts:
-                mat[ells] += numpy.dot(dwts[alpha], dexpansion_values[alpha].T)
+            dexpansion_values = expansion.tabulate(pts, max_deriv_order)
+            for Q in Qs_to_ells:
+                ells = Qs_to_ells[Q]
+                cur_pts = Qs_to_pts[Q]
+                indices = list(map(pts.index, cur_pts))
+                wshape = (len(ells), *tshape, len(cur_pts))
+                dwts = {alpha: numpy.zeros(wshape, "d") for alpha in dexpansion_values if sum(alpha) > 0}
+                if Q is None:
+                    for i, k in enumerate(ells):
+                        ell = self.nodes[k]
+                        for pt, wac_list in ell.deriv_dict.items():
+                            j = cur_pts.index(pt)
+                            for (w, alpha, c) in wac_list:
+                                dwts[alpha][i][c][j] = w
+                else:
+                    for i, k in enumerate(ells):
+                        ell = self.nodes[k]
+                        for alpha in ell.weights:
+                            dwts[alpha][i][ell.comp][:] = ell.weights[alpha]
+                for alpha in dwts:
+                    wts = dwts[alpha]
+                    expansion_values = dexpansion_values[alpha].T
+                    mat[ells] += numpy.dot(wts, expansion_values[indices])
         return mat
 
     def get_indices(self, restriction_domain, take_closure=True):