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.c2_elements import AlfeldC2, BrambleZlamalC2
 from FIAT.alfeld_sorokina import AlfeldSorokina
 from FIAT.arnold_qin import ArnoldQin
 from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div
@@ -94,6 +95,8 @@
                       "Hsieh-Clough-Tocher": HsiehCloughTocher,
                       "QuadraticPowellSabin6": QuadraticPowellSabin6,
                       "QuadraticPowellSabin12": QuadraticPowellSabin12,
+                      "Alfeld C2": AlfeldC2,
+                      "Bramble-Zlamal C2": BrambleZlamalC2,
                       "Alfeld-Sorokina": AlfeldSorokina,
                       "Arnold-Qin": ArnoldQin,
                       "Christiansen-Hu": ChristiansenHu,

FIAT/c2_elements.py

@@ -0,0 +1,116 @@
+# 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 C2DualSet(dual_set.DualSet):
+    """A DualSet for triangular C2 elements.
+
+    By default, this imposes C4 continuity at vertices for non-macroelements,
+    and the minimal C2 continuity for macroelements.
+    """
+    def __init__(self, ref_complex, degree, vorder=None, reduced=False, quad_scheme=None):
+        if vorder is None:
+            vorder = 2 if ref_complex.is_macrocell() else 4
+
+        if degree < 2*vorder+1:
+            raise ValueError(f"{type(self).__name__} only defined for degree >= {2*vorder+1}")
+
+        ref_el = ref_complex.get_parent() or ref_complex
+        if ref_el.get_shape() != TRIANGLE:
+            raise ValueError(f"{type(self).__name__} 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)}
+
+        # vorder jet at vertices
+        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, vorder+1) for alpha in mis(sd, i))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        k = degree - 2*vorder
+        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:
+            # Integral moments of normal derivatives against Jacobi polynomials along edges
+            phis = eval_jacobi_batch(vorder, vorder, k, xref)
+            dphis = 2*eval_jacobi_deriv_batch(vorder, vorder, k, xref, order=1)
+            ddphis = 4*eval_jacobi_deriv_batch(vorder, vorder, 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)))
+
+            # Interior moments against a basis for Pq
+            q = degree - 3 * (vorder // 2 + 1)
+            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 BrambleZlamalC2(finite_element.CiarletElement):
+    """The Bramble-Zlamal C2 element."""
+    def __init__(self, ref_el, degree=9, reduced=False, quad_scheme=None):
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
+        dual = C2DualSet(ref_el, degree, reduced=reduced, quad_scheme=quad_scheme)
+        super().__init__(poly_set, dual, degree, formdegree=0)
+
+
+def AlfeldC2Space(ref_el, degree):
+    """Construct the generalization of the quintic C2 spline on the double Alfeld split."""
+    # Construct the double Alfeld split by splitting twice
+    ref_complex = macro.AlfeldSplit(macro.AlfeldSplit(ref_el))
+    # 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
+    return macro.CkPolynomialSet(ref_complex, degree, order=order, variant="bubble")
+
+
+class AlfeldC2(finite_element.CiarletElement):
+    """The Alfeld C^2 macroelement on a double barycentric split.
+    See Section 7.5 of Lai & Schumacher for the quintic C^2 spline.
+    """
+    def __init__(self, ref_el, degree=5, reduced=False, quad_scheme=None):
+        poly_set = AlfeldC2Space(ref_el, degree)
+        ref_complex = poly_set.get_reference_element()
+        dual = C2DualSet(ref_complex, degree, reduced=reduced, quad_scheme=quad_scheme)
+        super().__init__(poly_set, dual, degree, formdegree=0)

FIAT/jacobi.py

@@ -82,17 +82,21 @@ def eval_jacobi_deriv(a, b, n, x):
         return 0.5 * (a + b + n + 1) * eval_jacobi(a + 1, b + 1, n - 1, x)
 
 
-def eval_jacobi_deriv_batch(a, b, n, xs):
-    """Evaluates the first derivatives of all jacobi polynomials with
+def eval_jacobi_deriv_batch(a, b, n, xs, order=1):
+    """Evaluates the derivative of the given order of all jacobi polynomials with
     weights a,b up to degree n.  xs is a numpy.array of points.
     Returns a two-dimensional array of points, where the
     rows correspond to the Jacobi polynomials and the
     columns correspond to the points."""
     results = numpy.zeros((n + 1, len(xs)), xs.dtype)
-    if n == 0:
+    if n+1 <= order:
         return results
     else:
-        results[1:, :] = eval_jacobi_batch(a + 1, b + 1, n - 1, xs)
-    for j in range(1, n + 1):
-        results[j, :] *= 0.5 * (a + b + j + 1)
+        results[order:, :] = eval_jacobi_batch(a + order, b + order, n - order, xs)
+    for j in range(order, n + 1):
+        z = 1
+        f = a + b + j + 1
+        for l in range(order):
+            z *= 0.5 * (f + l)
+        results[j, :] *= z
     return results

finat/__init__.py

@@ -23,6 +23,7 @@
 from .hct import HsiehCloughTocher, ReducedHsiehCloughTocher       # noqa: F401
 from .arnold_qin import ArnoldQin, ReducedArnoldQin                # noqa: F401
 from .christiansen_hu import ChristiansenHu                        # noqa: F401
+from .c2_elements import AlfeldC2, BrambleZlamalC2                 # noqa: F401
 from .alfeld_sorokina import AlfeldSorokina                        # noqa: F401
 from .guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanBubble, GuzmanNeilanH1div  # noqa: F401
 from .powell_sabin import QuadraticPowellSabin6, QuadraticPowellSabin12  # noqa: F401

finat/argyris.py

@@ -1,44 +1,66 @@
 import numpy
 from math import comb
+from itertools import chain
 
 import FIAT
 
-from gem import Literal, ListTensor
+from gem import Literal, ListTensor, Zero
 
 from finat.citations import cite
 from finat.fiat_elements import ScalarFiatElement
 from finat.physically_mapped import identity, PhysicallyMappedElement
 
 
+def _jet_transform(J, order):
+    """Basis transformation for derivative evaluation."""
+    if order == 0:
+        return identity(1)
+    sd = J.shape[0]
+    shape = (sd,)*order
+
+    # Mapping from multiindices to linearly-independent (flattened) components
+    mapping = {}
+    alphas = []
+    for indices in numpy.ndindex(shape):
+        alpha = [0] * sd
+        for i in indices:
+            alpha[i] += 1
+        alpha = tuple(alpha)
+        if alpha not in alphas:
+            alphas.append(alpha)
+        mapping[indices] = alphas.index(alpha)
+    # Inverse mapping
+    imapping = {v: k for k, v in mapping.items()}
+
+    # Get the transformation for a covariant tensor.
+    # We take the outer product, as each index maps with the Jacobian.
+    Jnp = numpy.asarray([[J[i, j] for j in range(sd)] for i in range(sd)])
+    Jprod = Jnp
+    for i in range(1, order):
+        Jprod = Jprod[..., None, None] * Jnp
+
+    # Deal with symmetries by contracting along linearly-dependent components.
+    B = numpy.full((len(alphas), len(alphas)), Zero(), dtype=object)
+    for i, ii in imapping.items():
+        for jj, j in mapping.items():
+            B[i, j] += Jprod[tuple(chain.from_iterable(zip(jj, ii)))]
+    return B
+
+
 def _vertex_transform(V, vorder, fiat_cell, coordinate_mapping):
-    """Basis transformation for evaluation, gradient, and hessian at vertices."""
+    """Basis transformation for jet at vertices."""
     sd = fiat_cell.get_spatial_dimension()
     top = fiat_cell.get_topology()
     bary, = fiat_cell.make_points(sd, 0, sd+1)
     J = coordinate_mapping.jacobian_at(bary)
 
-    gdofs = sd
-    G = [[J[j, i] for j in range(sd)] for i in range(sd)]
-
-    if vorder < 2:
-        hdofs = 0
-        H = [[]]
-    else:
-        hdofs = (sd*(sd+1))//2
-        indices = [(i, j) for i in range(sd) for j in range(i, sd)]
-        H = numpy.zeros((hdofs, hdofs), dtype=object)
-        for p, (i, j) in enumerate(indices):
-            for q, (m, n) in enumerate(indices):
-                H[p, q] = J[m, i] * J[n, j] + J[m, j] * J[n, i]
-        H[:, [i == j for i, j in indices]] *= 0.5
-
+    jet = [_jet_transform(J, k) for k in range(vorder+1)]
     s = 0
     for v in sorted(top[0]):
-        s += 1
-        V[s:s+gdofs, s:s+gdofs] = G
-        s += gdofs
-        V[s:s+hdofs, s:s+hdofs] = H
-        s += hdofs
+        for B in jet:
+            ndofs = len(B)
+            V[s:s+ndofs, s:s+ndofs] = B
+            s += ndofs
     return V
 
 

finat/c2_elements.py

@@ -0,0 +1,126 @@
+from finat.fiat_elements import ScalarFiatElement
+from finat.physically_mapped import identity, PhysicallyMappedElement
+from finat.citations import cite
+from finat.argyris import (_jet_transform, _vertex_transform,
+                           _normal_tangential_transform)
+from gem import ListTensor
+
+import FIAT
+import numpy
+from math import comb
+
+
+class C2Element(PhysicallyMappedElement):
+
+    def basis_transformation(self, coordinate_mapping):
+        top = self.cell.topology
+        sd = self.cell.get_spatial_dimension()
+        entity_ids = self._element.entity_dofs()
+
+        # Detect the maximum derivative order of the vertex dofs
+        nodes = self._element.dual_basis()
+        vorder = max(nodes[i].max_deriv_order for i in entity_ids[0][0])
+
+        V = identity(self.space_dimension())
+        _vertex_transform(V, vorder, self.cell, coordinate_mapping)
+
+        bary, = self.cell.make_points(sd, 0, sd+1)
+        J = coordinate_mapping.jacobian_at(bary)
+        detJ = coordinate_mapping.detJ_at(bary)
+        Thetainv = _jet_transform(J, 2)
+
+        ns = coordinate_mapping.physical_normals()
+        ts = coordinate_mapping.physical_tangents()
+        lens = coordinate_mapping.physical_edge_lengths()
+        nhats = coordinate_mapping.reference_normals()
+        thats = coordinate_mapping.normalized_reference_edge_tangents()
+
+        n0 = self.degree - 2*vorder - 1
+        n1 = n0 + 1
+        for e in top[1]:
+            v0, v1 = top[1][e]
+            vid0 = entity_ids[0][v0]
+            vid1 = entity_ids[0][v1]
+            eids = entity_ids[1][e]
+            emoments = (eids[:n0], eids[n0:n0+n1], eids[n0+n1:])
+
+            G = numpy.array([[u[e, j] for j in range(sd)] for u in (ns, ts)])
+            Ghat = numpy.array([[u[e, j] for j in range(sd)] for u in (nhats, thats)])
+            Gamma = _jet_transform(G, 2)
+            Gammainvhat = _jet_transform(Ghat.T, 2)
+
+            B2 = (Gammainvhat @ Thetainv) @ Gamma
+            beta = B2[0, 1:] @ G / lens[e]
+
+            Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e)
+            if self.avg:
+                Bnn = Bnn * lens[e]
+
+            # first derivative moments
+            for k, s1 in enumerate(emoments[1], start=1):
+                # Derivative of Jacobi polynomial at the endpoints
+                dP1 = comb(k + vorder, k-1) * (2*vorder+k+1)
+                dP0 = (-1)**k * dP1
+
+                V[s1, s1] = Bnn
+                V[s1, vid0[0]] = dP0 * Bnt
+                V[s1, vid1[0]] = dP1 * Bnt
+                if k > 1:
+                    s0 = emoments[0][k-2]
+                    V[s1, s0] = -1 * Bnt
+
+            # second derivative moments
+            for k, s2 in enumerate(emoments[2]):
+                # Jacobi polynomial at the endpoints
+                P1 = comb(k + vorder, k)
+                P0 = -(-1)**k * P1
+
+                V[s2, s2] = B2[0, 0]
+                V[s2, vid0[1:sd+1]] = P0 * beta
+                V[s2, vid1[1:sd+1]] = P1 * beta
+                if k > 0:
+                    s1 = emoments[1][k-1]
+                    V[s2, s1] = -2 * Bnt * V[s1, s1]
+                    V[s2, vid0[0]] = -1 * Bnt * V[s1, vid0[0]]
+                    V[s2, vid1[0]] = -1 * Bnt * V[s1, vid1[0]]
+                if k > 1:
+                    s0 = emoments[0][k-2]
+                    V[s2, s0] = -1 * Bnt * V[s1, s0]
+
+        # Now let's fix the scaling.
+        h = coordinate_mapping.cell_size()
+        for v in top[0]:
+            vids = entity_ids[0][v]
+            scale = 1 / h[v]
+            F = scale
+            iend = 1
+            # scale the jet at vertices
+            for k in range(1, vorder+1):
+                istart = iend
+                iend = istart + comb(k+sd-1, sd-1)
+                V[:, vids[istart:iend]] *= F
+                F *= scale
+
+        for e in top[1]:
+            eids = entity_ids[1][e]
+            emoments = (eids[:n0], eids[n0:n0+n1], eids[n0+n1:])
+            he = (1/len(top[1][e])) * sum(h[v] for v in top[1][e])
+            # scale first and second derivative moments
+            V[:, emoments[1]] *= 1 / he
+            V[:, emoments[2]] *= 1 / (he * he)
+
+        return ListTensor(V.T)
+
+
+class BrambleZlamalC2(C2Element, ScalarFiatElement):
+    def __init__(self, cell, degree=9, avg=True):
+        cite("BrambleZlamal1970")
+        self.avg = avg
+        super().__init__(FIAT.BrambleZlamalC2(cell, degree))
+
+
+class AlfeldC2(C2Element, ScalarFiatElement):
+    def __init__(self, cell, degree=5, avg=True):
+        cite("Alfeld1984")
+        self.avg = avg
+        super().__init__(FIAT.AlfeldC2(cell, degree))