Walkington C1 macroelement

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.walkington import Walkington
 from FIAT.histopolation import Histopolation
 from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite  # noqa: F401
 
@@ -117,7 +118,8 @@
                       "Conforming Arnold-Winther": ArnoldWinther,
                       "Nonconforming Arnold-Winther": ArnoldWintherNC,
                       "Hu-Zhang": HuZhang,
-                      "Mardal-Tai-Winther": MardalTaiWinther}
+                      "Mardal-Tai-Winther": MardalTaiWinther,
+                      "Walkington": Walkington}
 
 # List of extra elements
 extra_elements = {"P0": P0}

FIAT/bell.py

@@ -9,68 +9,51 @@
 # bfs, but the extra three are used in the transformation theory.
 
 from FIAT import finite_element, polynomial_set, dual_set, functional
-from FIAT.reference_element import TRIANGLE, ufc_simplex
+from FIAT.reference_element import TRIANGLE
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.jacobi import eval_jacobi
 
 
 class BellDualSet(dual_set.DualSet):
-    def __init__(self, ref_el):
-        entity_ids = {}
-        nodes = []
-        cur = 0
-
-        # make nodes by getting points
-        # need to do this dimension-by-dimension, facet-by-facet
+    def __init__(self, ref_el, degree):
         top = ref_el.get_topology()
-        verts = ref_el.get_vertices()
         sd = ref_el.get_spatial_dimension()
-        if ref_el.get_shape() != TRIANGLE:
-            raise ValueError("Bell only defined on triangles")
-
-        pd = functional.PointDerivative
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
         # get jet at each vertex
-
-        entity_ids[0] = {}
         for v in sorted(top[0]):
-            nodes.append(functional.PointEvaluation(ref_el, verts[v]))
-
-            # first derivatives
-            for i in range(sd):
-                alpha = [0] * sd
-                alpha[i] = 1
-                nodes.append(pd(ref_el, verts[v], alpha))
-
-            # second derivatives
-            alphas = [[2, 0], [1, 1], [0, 2]]
-            for alpha in alphas:
-                nodes.append(pd(ref_el, verts[v], alpha))
+            cur = len(nodes)
+            x, = ref_el.make_points(0, v, degree)
+            nodes.append(functional.PointEvaluation(ref_el, x))
 
-            entity_ids[0][v] = list(range(cur, cur + 6))
-            cur += 6
+            # first and second derivatives
+            nodes.extend(functional.PointDerivative(ref_el, x, alpha)
+                         for i in (1, 2) for alpha in polynomial_set.mis(sd, i))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
 
         # we need an edge quadrature rule for the moment
-        from FIAT.quadrature_schemes import create_quadrature
-        from FIAT.jacobi import eval_jacobi
-        rline = ufc_simplex(1)
-        q1d = create_quadrature(rline, 8)
-        q1dpts = q1d.get_points()[:, 0]
-        leg4_at_qpts = eval_jacobi(0, 0, 4, 2.0*q1dpts - 1)
+        facet = ref_el.construct_subelement(1)
+        Q_ref = create_quadrature(facet, 2*(degree-1))
+        x = facet.compute_barycentric_coordinates(Q_ref.get_points())
+        leg4_at_qpts = eval_jacobi(0, 0, 4, x[:, 1] - x[:, 0])
 
-        imond = functional.IntegralMomentOfNormalDerivative
-        entity_ids[1] = {}
         for e in sorted(top[1]):
-            entity_ids[1][e] = [18+e]
-            nodes.append(imond(ref_el, e, q1d, leg4_at_qpts))
-
-        entity_ids[2] = {0: []}
+            cur = len(nodes)
+            nodes.append(functional.IntegralMomentOfNormalDerivative(ref_el, e, Q_ref, leg4_at_qpts))
+            entity_ids[1][e].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
 
 class Bell(finite_element.CiarletElement):
     """The Bell finite element."""
 
-    def __init__(self, ref_el):
-        poly_set = polynomial_set.ONPolynomialSet(ref_el, 5)
-        dual = BellDualSet(ref_el)
-        super().__init__(poly_set, dual, 5)
+    def __init__(self, ref_el, degree=5):
+        if ref_el.get_shape() != TRIANGLE:
+            raise ValueError(f"{type(self).__name__} only defined on triangles")
+        if degree != 5:
+            raise ValueError(f"{type(self).__name__} only defined for degree = 5.")
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
+        dual = BellDualSet(ref_el, degree)
+        super().__init__(poly_set, dual, degree)

FIAT/functional.py

@@ -342,8 +342,9 @@ def __init__(self, ref_el, s, Q, f_at_qpts, comp=(), shp=()):
 class IntegralMomentOfNormalDerivative(Functional):
     """Functional giving normal derivative integrated against some function on a facet."""
 
-    def __init__(self, ref_el, facet_no, Q, f_at_qpts):
-        n = ref_el.compute_normal(facet_no)
+    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
@@ -364,6 +365,39 @@ def __init__(self, ref_el, facet_no, Q, f_at_qpts):
                          {}, 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")
+
+
 class FrobeniusIntegralMoment(IntegralMoment):
 
     def __init__(self, ref_el, Q, f_at_qpts, nm=None):

FIAT/walkington.py

@@ -0,0 +1,120 @@
+# Copyright (C) 2025 Pablo D. Brubeck
+#
+# This file is part of FIAT (https://www.fenicsproject.org)
+#
+# SPDX-License-Identifier:    LGPL-3.0-or-later
+
+# This is not quite Walkington, but is 65-dofs and includes 20 extra constraint
+# functionals.  The first 45 basis functions are the reference element
+# bfs, but the extra 20 are used in the transformation theory.
+
+from FIAT import finite_element, polynomial_set, dual_set, macro
+from FIAT.expansions import polynomial_dimension
+from FIAT.functional import (
+    PointEvaluation, PointDerivative,
+    IntegralMomentOfNormalDerivative,
+    IntegralMomentOfBidirectionalDerivative,
+)
+from FIAT.reference_element import TETRAHEDRON
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.jacobi import eval_jacobi
+import numpy
+
+
+class WalkingtonDualSet(dual_set.DualSet):
+    def __init__(self, ref_el, degree):
+        top = ref_el.get_topology()
+        sd = ref_el.get_spatial_dimension()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
+
+        # Vertex dofs: second order jet
+        for v in sorted(top[0]):
+            cur = len(nodes)
+            x, = ref_el.make_points(0, v, degree)
+            nodes.append(PointEvaluation(ref_el, x))
+
+            # first and second derivatives
+            nodes.extend(PointDerivative(ref_el, x, alpha)
+                         for i in (1, 2) for alpha in polynomial_set.mis(sd, i))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        # Face dofs: moments or normal derivative
+        ref_face = ref_el.construct_subelement(2)
+        Q_face = create_quadrature(ref_face, degree-1)
+        f_at_qpts = numpy.ones(Q_face.get_weights().shape)
+        for face in sorted(top[2]):
+            cur = len(nodes)
+            nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, f_at_qpts))
+            entity_ids[2][face].extend(range(cur, len(nodes)))
+
+        # Interior dof: point evaluation at barycenter
+        for entity in top[sd]:
+            cur = len(nodes)
+            x, = ref_el.make_points(sd, entity, sd+1)
+            nodes.append(PointEvaluation(ref_el, x))
+            entity_ids[sd][entity].extend(range(cur, len(nodes)))
+
+        # Constraint dofs
+        # Face-edge constraint: normal derivative along edge is cubic
+        edges = ref_el.get_connectivity()[(2, 1)]
+        ref_edge = ref_el.construct_subelement(1)
+        Q_edge = create_quadrature(ref_edge, 2*(degree-1))
+        x = ref_edge.compute_barycentric_coordinates(Q_edge.get_points())
+        leg4_at_qpts = eval_jacobi(0, 0, 4, x[:, 1] - x[:, 0])
+        # Face constraint: normal derivative is cubic
+        Q_face, phi = face_constraint(ref_face)
+
+        for face in sorted(top[2]):
+            cur = len(nodes)
+            thats = ref_el.compute_tangents(sd-1, face)
+            nface = -numpy.cross(*thats)
+            nface /= numpy.linalg.norm(nface)
+
+            for i, e in enumerate(edges[face]):
+                Q = FacetQuadratureRule(ref_face, 1, i, Q_edge)
+
+                te = ref_el.compute_edge_tangent(e)
+                nfe = numpy.cross(te, nface)
+                nfe /= numpy.linalg.norm(nfe)
+                nfe /= Q.jacobian_determinant()
+                nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q, leg4_at_qpts, n=nfe))
+
+            Q = FacetQuadratureRule(ref_el, 2, face, Q_face)
+            nface /= Q.jacobian_determinant()
+            nodes.extend(IntegralMomentOfBidirectionalDerivative(ref_el, Q, phi, nface, t) for t in thats)
+            entity_ids[2][face].extend(range(cur, len(nodes)))
+
+        super().__init__(nodes, ref_el, entity_ids)
+
+
+class Walkington(finite_element.CiarletElement):
+    """The Walkington C1 macroelement."""
+
+    def __init__(self, ref_el, degree=5):
+        if ref_el.get_shape() != TETRAHEDRON:
+            raise ValueError(f"{type(self).__name__} only defined on tetrahedron")
+        if degree != 5:
+            raise ValueError(f"{type(self).__name__} only defined for degree=5.")
+
+        dual = WalkingtonDualSet(ref_el, degree)
+        ref_complex = macro.AlfeldSplit(ref_el)
+        poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1, vorder=4, variant="bubble")
+        super().__init__(poly_set, dual, degree)
+
+
+def face_constraint(ref_face):
+    k = 3
+    sd = ref_face.get_spatial_dimension()
+    Q_face = create_quadrature(ref_face, 2*k)
+    dimPkm1 = polynomial_dimension(ref_face, k-1)
+
+    pts = list(Q_face.get_points()[:3])
+    pts.append(Q_face.get_points()[-1])
+    P = polynomial_set.ONPolynomialSet(ref_face, k)
+    Pk = P.tabulate(pts)[(0,)*sd][dimPkm1:]
+    c = numpy.linalg.solve(Pk.T, [0, 0, 0, 1])
+    Pk = P.tabulate(Q_face.get_points())[(0,)*sd][dimPkm1:]
+    phi = numpy.dot(c, Pk)
+    return Q_face, phi

finat/__init__.py

@@ -29,6 +29,7 @@
 from .johnson_mercier import JohnsonMercier                        # noqa: F401
 from .mtw import MardalTaiWinther                                  # noqa: F401
 from .morley import Morley                                         # noqa: F401
+from .walkington import Walkington                                 # noqa: F401
 from .direct_serendipity import DirectSerendipity                  # noqa: F401
 from .tensorfiniteelement import TensorFiniteElement               # noqa: F401
 from .tensor_product import TensorProductElement                   # noqa: F401

finat/argyris.py

@@ -42,10 +42,16 @@ def _vertex_transform(V, vorder, fiat_cell, coordinate_mapping):
     return V
 
 
-def _normal_tangential_transform(fiat_cell, J, detJ, f):
-    R = numpy.array([[0, 1], [-1, 0]])
-    that = fiat_cell.compute_edge_tangent(f)
-    nhat = R @ that
+def _normal_tangential_transform(fiat_cell, J, detJ, edge, face=None):
+    that = fiat_cell.compute_edge_tangent(edge)
+    if fiat_cell.get_spatial_dimension() == 2:
+        R = numpy.array([[0, 1], [-1, 0]])
+        nhat = R @ that
+    else:
+        nface = fiat_cell.compute_scaled_normal(face)
+        nface /= numpy.linalg.norm(nface)
+        nhat = numpy.cross(that, nface)
+
     Jn = J @ Literal(nhat)
     Jt = J @ Literal(that)
     alpha = Jn @ Jt

finat/bell.py

@@ -12,7 +12,7 @@
 class Bell(PhysicallyMappedElement, ScalarFiatElement):
     def __init__(self, cell, degree=5):
         cite("Bell1969")
-        super().__init__(FIAT.Bell(cell))
+        super().__init__(FIAT.Bell(cell, degree=degree))
 
         reduced_dofs = deepcopy(self._element.entity_dofs())
         sd = cell.get_spatial_dimension()

finat/citations.py

@@ -259,3 +259,13 @@ def cite(*args, **kwargs):
  doi={https://doi.org/10.1023/A:1004420829610},
 }
 """)
+    petsctools.add_citation("Walkington2010", """
+author = {Walkington, Noel J.},
+title = {{A \\$C^1\\$ Tetrahedral Finite Element without Edge Degrees of Freedom}},
+journal = {SIAM Journal on Numerical Analysis},
+volume = {52},
+number = {1},
+pages = {330-342},
+year = {2014},
+doi = {10.1137/130912013},
+}""")

finat/element_factory.py

@@ -90,6 +90,7 @@
                       "Nonconforming Arnold-Winther": finat.ArnoldWintherNC,
                       "Hu-Zhang": finat.HuZhang,
                       "Mardal-Tai-Winther": finat.MardalTaiWinther,
+                      "Walkington": finat.Walkington,
                       # These require special treatment
                       "Q": None,
                       "DQ": None,