Automate zany basis transformations

finat/__init__.py

@@ -43,6 +43,7 @@
 from .nodal_enriched import NodalEnrichedElement                   # noqa: F401
 from .quadrature_element import QuadratureElement, make_quadrature_element  # noqa: F401
 from .restricted import RestrictedElement                          # noqa: F401
+from .functional import Functional                                 # noqa: F401
 from .runtime_tabulated import RuntimeTabulated                    # noqa: F401
 from . import quadrature                                           # noqa: F401
 from . import cell_tools                                           # noqa: F401

finat/argyris.py

@@ -8,7 +8,9 @@
 
 from finat.citations import cite
 from finat.fiat_elements import ScalarFiatElement
-from finat.physically_mapped import identity, PhysicallyMappedElement
+from finat.physically_mapped import (identity, PhysicalGeometry,
+                                     PhysicallyMappedElement)
+from finat.zany import zany_basis_transformation
 
 
 def _jet_transform(J, order):
@@ -128,70 +130,21 @@ def _edge_transform(V, vorder, eorder, fiat_cell, coordinate_mapping, avg=False)
 
 
 class Argyris(PhysicallyMappedElement, ScalarFiatElement):
+    """The Argyris element.
+
+    The basis transformation is derived automatically from the FIAT
+    dual basis by :func:`finat.zany.zany_basis_transformation`.
+    """
     def __init__(self, cell, degree=5, variant=None, avg=False):
         cite("Argyris1968")
         if variant is None:
             variant = "integral"
         if variant == "point" and degree != 5:
             raise NotImplementedError("Degree must be 5 for 'point' variant of Argyris")
-        fiat_element = FIAT.Argyris(cell, degree, variant=variant)
         self.variant = variant
         self.avg = avg
-        super().__init__(fiat_element)
-
-    def basis_transformation(self, coordinate_mapping):
-        sd = self.cell.get_spatial_dimension()
-        top = self.cell.get_topology()
-
-        V = identity(self.space_dimension())
-
-        vorder = 2
-        voffset = comb(sd + vorder, vorder)
-        eorder = self.degree - 5
-
-        _vertex_transform(V, vorder, self.cell, coordinate_mapping)
-        if self.variant == "integral":
-            _edge_transform(V, vorder, eorder, self.cell, coordinate_mapping, avg=self.avg)
-        else:
-            bary, = self.cell.make_points(sd, 0, sd+1)
-            J = coordinate_mapping.jacobian_at(bary)
-            detJ = coordinate_mapping.detJ_at(bary)
-            pel = coordinate_mapping.physical_edge_lengths()
-            for e in sorted(top[1]):
-                s = len(top[0]) * voffset + e * (eorder+1)
-                v0id, v1id = (v * voffset for v in top[1][e])
-                Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e)
-
-                # edge midpoint normal derivative
-                V[s, s] = Bnn * pel[e]
-
-                # vertex points
-                V[s, v1id] = 15/8 * Bnt
-                V[s, v0id] = -V[s, v1id]
-
-                # vertex derivatives
-                for i in range(sd):
-                    V[s, v1id+1+i] = -7/16 * Bnt * Jt[i]
-                    V[s, v0id+1+i] = V[s, v1id+1+i]
-
-                # second derivatives
-                tau = [Jt[0]*Jt[0], 2*Jt[0]*Jt[1], Jt[1]*Jt[1]]
-                for i in range(len(tau)):
-                    V[s, v1id+3+i] = 1/32 * Bnt * tau[i]
-                    V[s, v0id+3+i] = -V[s, v1id+3+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])
-
-        if self.variant == "point":
-            eoffset = 2 * eorder + 1
-            for e in sorted(top[1]):
-                v0, v1 = top[1][e]
-                s = len(top[0]) * voffset + e * eoffset
-                V[:, s:s+eorder+1] *= 2 / (h[v0] + h[v1])
-
-        return ListTensor(V.T)
+        super().__init__(FIAT.Argyris(cell, degree, variant=variant))
+
+    def basis_transformation(self, coordinate_mapping: PhysicalGeometry) -> ListTensor:
+        return zany_basis_transformation(self._element, coordinate_mapping,
+                                         avg=self.avg)

finat/bell.py

@@ -1,15 +1,24 @@
 import FIAT
-from math import comb
+from copy import deepcopy
+
 from gem import ListTensor
 
 from finat.citations import cite
 from finat.fiat_elements import ScalarFiatElement
-from finat.physically_mapped import identity, PhysicallyMappedElement
-from finat.argyris import _vertex_transform, _normal_tangential_transform
-from copy import deepcopy
+from finat.physically_mapped import PhysicalGeometry, PhysicallyMappedElement
+from finat.zany import zany_basis_transformation
 
 
 class Bell(PhysicallyMappedElement, ScalarFiatElement):
+    """The Bell element.
+
+    FIAT provides the extended element on the full quintic space, with
+    the cubic normal derivative constraints appended as extra edge
+    nodes; the transformation of the extended element is derived
+    automatically by :func:`finat.zany.zany_basis_transformation`, and
+    the constraint degrees of freedom are dropped from the physical
+    element.
+    """
     def __init__(self, cell, degree=5):
         cite("Bell1969")
         super().__init__(FIAT.Bell(cell, degree=degree))
@@ -20,56 +29,12 @@ def __init__(self, cell, degree=5):
             reduced_dofs[sd-1][entity] = []
         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)
-
-        voffset = comb(sd + vorder, vorder)
-        for e in sorted(top[1]):
-            s = len(top[0]) * voffset + e
-            v0id, v1id = (v * voffset for v in top[1][e])
-            Bnn, Bnt, Jt = _normal_tangential_transform(self.cell, J, detJ, e)
-
-            # vertex points
-            V[s, v1id] = 1/21 * Bnt
-            V[s, v0id] = -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
-            tau = [Jt[0]*Jt[0], 2*Jt[0]*Jt[1], Jt[1]*Jt[1]]
-            for i in range(len(tau)):
-                V[s, v1id+3+i] = 1/252 * Bnt * tau[i]
-                V[s, v0id+3+i] = -V[s, v1id+3+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)
+    def basis_transformation(self, coordinate_mapping: PhysicalGeometry) -> ListTensor:
+        return zany_basis_transformation(self._element, coordinate_mapping,
+                                         ndof=self.space_dimension())
 
-    # This wipes out the edge dofs.  FIAT gives a 21 DOF element
-    # because we need some extra functions to help with transforming
-    # under the edge constraint.  However, we only have an 18 DOF
-    # element.
+    # The extended FIAT element has 21 basis functions, but the Bell
+    # element only keeps the 18 vertex degrees of freedom.
     def entity_dofs(self):
         return self._entity_dofs