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

@@ -349,8 +349,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

FIAT/walkington.py

@@ -0,0 +1,100 @@
+# 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.functional import PointEvaluation, PointDerivative, IntegralMomentOfNormalDerivative
+from FIAT.reference_element import TETRAHEDRON
+from FIAT.quadrature import FacetQuadratureRule
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.jacobi import eval_jacobi
+from FIAT.hierarchical import make_dual_bubbles
+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 against cubic bubble
+        ref_face = ref_el.construct_subelement(2)
+        Q_face, phis = make_dual_bubbles(ref_face, degree-1)
+
+        for face in sorted(top[2]):
+            cur = len(nodes)
+            nface = ref_el.compute_scaled_normal(face)
+            nface /= numpy.dot(nface, nface)
+
+            nodes.append(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phis[0], n=nface))
+            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
+        # Face constraint: normal derivative 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])
+
+        for face in sorted(top[2]):
+            cur = len(nodes)
+            nface = ref_el.compute_scaled_normal(face)
+            nface /= numpy.dot(nface, 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))
+
+            nodes.extend(IntegralMomentOfNormalDerivative(ref_el, face, Q_face, phi, n=nface) for phi in phis[1:])
+            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)

finat/__init__.py

@@ -25,6 +25,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 .spectral import GaussLobattoLegendre, GaussLegendre, Legendre, IntegratedLegendre, FDMLagrange, FDMQuadrature, FDMDiscontinuousLagrange, FDMBrokenH1, FDMBrokenL2, FDMHermite  # noqa: F401
 from .tensorfiniteelement import TensorFiniteElement               # noqa: F401

finat/argyris.py

@@ -41,10 +41,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):
         if Citations is not None:
             Citations().register("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/walkington.py

@@ -0,0 +1,90 @@
+import FIAT
+from FIAT.polynomial_set import mis
+from math import comb
+
+from gem import ListTensor
+
+from finat.fiat_elements import ScalarFiatElement
+from finat.physically_mapped import Citations, identity, PhysicallyMappedElement
+from finat.argyris import _vertex_transform, _normal_tangential_transform
+from copy import deepcopy
+
+
+class Walkington(PhysicallyMappedElement, ScalarFiatElement):
+    def __init__(self, cell, degree=5):
+        if Citations is not None:
+            Citations().register("Bell1969")
+        super().__init__(FIAT.Walkington(cell, degree=degree))
+
+        reduced_dofs = deepcopy(self._element.entity_dofs())
+        sd = cell.get_spatial_dimension()
+        for entity in reduced_dofs[sd-1]:
+            reduced_dofs[sd-1][entity] = reduced_dofs[sd-1][entity][:1]
+        self._entity_dofs = reduced_dofs
+
+    def basis_transformation(self, coordinate_mapping):
+        # Jacobian at barycenter
+        sd = self.cell.get_spatial_dimension()
+        top = self.cell.get_topology()
+        bary, = self.cell.make_points(sd, 0, sd+1)
+        J = coordinate_mapping.jacobian_at(bary)
+        detJ = coordinate_mapping.detJ_at(bary)
+
+        numbf = self._element.space_dimension()
+        ndof = self.space_dimension()
+        # rectangular to toss out the constraint dofs
+        V = identity(numbf, ndof)
+
+        vorder = 2
+        _vertex_transform(V, vorder, self.cell, coordinate_mapping)
+
+        offset = ndof
+        voffset = comb(sd + vorder, vorder)
+        foffset = len(self._element.entity_dofs()[2][0]) - len(self.entity_dofs()[2][0])
+
+        edges = self.cell.get_connectivity()[(2, 1)]
+        for f in sorted(top[2]):
+            q = len(top[0]) * voffset + f
+            V[q, q] *= detJ
+
+            for j, e in enumerate(edges[f]):
+                s = offset + foffset * f + j
+
+                v0id, v1id = (v * voffset for v in top[1][e])
+                Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e, face=f)
+
+                # vertex points
+                V[s, v1id] = 1/21 * Bnt
+                V[s, v0id] = -1 * V[s, v1id]
+
+                # vertex derivatives
+                for i in range(sd):
+                    V[s, v1id+1+i] = -1/42 * Bnt * Jt[i]
+                    V[s, v0id+1+i] = V[s, v1id+1+i]
+
+                # second derivatives
+                for i, alpha in enumerate(mis(sd, 2)):
+                    ids = tuple(k for k, ak in enumerate(alpha) if ak)
+                    a, b = ids[0], ids[-1]
+                    tau = (1 + (a != b)) * Jt[a] * Jt[b]
+                    V[s, v1id+sd+1+i] = 1/252 * Bnt * tau
+                    V[s, v0id+sd+1+i] = -1 * V[s, v1id+sd+1+i]
+
+        # Patch up conditioning
+        h = coordinate_mapping.cell_size()
+        for v in sorted(top[0]):
+            s = voffset * v + 1
+            V[:, s:s+sd] *= 1/h[v]
+            V[:, s+sd:voffset*(v+1)] *= 1/(h[v]*h[v])
+
+        return ListTensor(V.T)
+
+    # This wipes out the edge dofs.  FIAT gives a 65 DOF element
+    # because we need some extra functions to help with transforming
+    # under the edge constraint.  However, we only have an 45 DOF
+    # element.
+    def entity_dofs(self):
+        return self._entity_dofs
+
+    def space_dimension(self):
+        return 45

test/finat/test_zany_mapping.py

@@ -104,6 +104,7 @@ def test_C1_triangle(ref_to_phys, element):
 
 @pytest.mark.parametrize("element", [
                          finat.Morley,
+                         finat.Walkington,
                          ])
 def test_C1_tetrahedron(ref_to_phys, element):
     check_zany_mapping(element, ref_to_phys[3])