C0/C1 Macroelements (#123)

* DG variant in constructor
* Support FIAT macro elements
* Legendre variant
* TensorFiniteElement overload complex
* cleanup Bell, register and attempt to transform HCT

---------

Co-authored-by: Rob Kirby <robert.c.kirby@gmail.com>

finat/__init__.py

@@ -11,6 +11,7 @@
 from .fiat_elements import FacetBubble  # noqa: F401
 from .fiat_elements import KongMulderVeldhuizen  # noqa: F401
 from .argyris import Argyris            # noqa: F401
+from .hct import HsiehCloughTocher, ReducedHsiehCloughTocher   # noqa: F401
 from .bell import Bell                  # noqa: F401
 from .hermite import Hermite            # noqa: F401
 from .mtw import MardalTaiWinther       # noqa: F401
@@ -32,3 +33,4 @@
 from .restricted import RestrictedElement          # noqa: F401
 from .runtime_tabulated import RuntimeTabulated  # noqa: F401
 from . import quadrature  # noqa: F401
+from . import cell_tools  # noqa: F401

finat/cell_tools.py

@@ -0,0 +1,5 @@
+"""Find the maximal complex in a list of cell complexes.
+This is a pass-through from FIAT so that FInAT clients
+(e.g. tsfc) don't have to directly import FIAT."""
+
+from FIAT.reference_element import max_complex  # noqa: F401

finat/cube.py

@@ -1,9 +1,9 @@
-from __future__ import absolute_import, print_function, division
+from __future__ import absolute_import, division, print_function
 
-from FIAT.reference_element import UFCHexahedron, UFCQuadrilateral
-from FIAT.reference_element import compute_unflattening_map, flatten_entities, flatten_permutations
+from FIAT.reference_element import (UFCHexahedron, UFCQuadrilateral,
+                                    compute_unflattening_map, flatten_entities,
+                                    flatten_permutations)
 from FIAT.tensor_product import FlattenedDimensions as FIAT_FlattenedDimensions
-
 from gem.utils import cached_property
 
 from finat.finiteelementbase import FiniteElementBase
@@ -29,6 +29,10 @@ def cell(self):
         else:
             raise NotImplementedError("Cannot guess cell for spatial dimension %s" % dim)
 
+    @property
+    def complex(self):
+        return self.product.complex
+
     @property
     def degree(self):
         unique_degree, = set(self.product.degree)

finat/direct_serendipity.py

@@ -33,6 +33,10 @@ def __init__(self, cell, degree):
     def cell(self):
         return self._cell
 
+    @property
+    def complex(self):
+        return self._cell
+
     @property
     def degree(self):
         return self._degree

finat/discontinuous.py

@@ -16,6 +16,10 @@ def __init__(self, element):
     def cell(self):
         return self.element.cell
 
+    @property
+    def complex(self):
+        return self.element.complex
+
     @property
     def degree(self):
         return self.element.degree

finat/enriched.py

@@ -1,12 +1,10 @@
 from functools import partial
-from operator import add, methodcaller
 from itertools import chain
-
-import numpy
+from operator import add, methodcaller
 
 import FIAT
-
 import gem
+import numpy
 from gem.utils import cached_property
 
 from finat.finiteelementbase import FiniteElementBase
@@ -30,6 +28,10 @@ def cell(self):
         result, = set(elem.cell for elem in self.elements)
         return result
 
+    @cached_property
+    def complex(self):
+        return FIAT.reference_element.max_complex(set(elem.complex for elem in self.elements))
+
     @cached_property
     def degree(self):
         return tree_map(max, *[elem.degree for elem in self.elements])

finat/fiat_elements.py

@@ -1,9 +1,7 @@
-import numpy as np
-import sympy as sp
-
 import FIAT
-
 import gem
+import numpy as np
+import sympy as sp
 from gem.utils import cached_property
 
 from finat.finiteelementbase import FiniteElementBase
@@ -53,6 +51,10 @@ def __init__(self, fiat_element):
     def cell(self):
         return self._element.get_reference_element()
 
+    @property
+    def complex(self):
+        return self._element.get_reference_complex()
+
     @property
     def degree(self):
         # Requires FIAT.CiarletElement
@@ -120,21 +122,21 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None):
 
             exprs = []
             for table in table_roll:
-                if derivative < self.degree:
+                if derivative == self.degree and not self.complex.is_macrocell():
+                    # Make sure numerics satisfies theory
+                    exprs.append(gem.Literal(table[0]))
+                elif derivative > self.degree:
+                    # Make sure numerics satisfies theory
+                    assert np.allclose(table, 0.0)
+                    exprs.append(gem.Literal(np.zeros(self.index_shape)))
+                else:
                     point_indices = ps.indices
                     point_shape = tuple(index.extent for index in point_indices)
 
                     exprs.append(gem.partial_indexed(
                         gem.Literal(table.reshape(point_shape + index_shape)),
                         point_indices
                     ))
-                elif derivative == self.degree:
-                    # Make sure numerics satisfies theory
-                    exprs.append(gem.Literal(table[0]))
-                else:
-                    # Make sure numerics satisfies theory
-                    assert np.allclose(table, 0.0)
-                    exprs.append(gem.Literal(np.zeros(self.index_shape)))
             if self.value_shape:
                 # As above, this extent may be different from that
                 # advertised by the finat element.
@@ -362,8 +364,8 @@ def __init__(self, cell, degree):
 
 
 class DiscontinuousLagrange(ScalarFiatElement):
-    def __init__(self, cell, degree):
-        super(DiscontinuousLagrange, self).__init__(FIAT.DiscontinuousLagrange(cell, degree))
+    def __init__(self, cell, degree, variant=None):
+        super(DiscontinuousLagrange, self).__init__(FIAT.DiscontinuousLagrange(cell, degree, variant=variant))
 
 
 class Real(DiscontinuousLagrange):

finat/finiteelementbase.py

@@ -1,9 +1,8 @@
-from abc import ABCMeta, abstractproperty, abstractmethod
+from abc import ABCMeta, abstractmethod, abstractproperty
 from itertools import chain
 
-import numpy
-
 import gem
+import numpy
 from gem.interpreter import evaluate
 from gem.optimise import delta_elimination, sum_factorise, traverse_product
 from gem.utils import cached_property
@@ -17,6 +16,11 @@ class FiniteElementBase(metaclass=ABCMeta):
     def cell(self):
         '''The reference cell on which the element is defined.'''
 
+    @property
+    def complex(self):
+        '''The reference cell complex over which bases are defined.
+           Can be different than self.cell in the case of macro elements.'''
+
     @abstractproperty
     def degree(self):
         '''The degree of the embedding polynomial space.

finat/hct.py

@@ -0,0 +1,147 @@
+import FIAT
+import numpy
+from gem import ListTensor, Literal, partial_indexed
+
+from finat.fiat_elements import ScalarFiatElement
+from finat.physically_mapped import Citations, PhysicallyMappedElement
+from copy import deepcopy
+
+
+class HsiehCloughTocher(PhysicallyMappedElement, ScalarFiatElement):
+    def __init__(self, cell, degree, avg=False):
+        if degree != 3:
+            raise ValueError("Degree must be 3 for HCT element")
+        if Citations is not None:
+            Citations().register("Clough1965")
+        self.avg = avg
+        super().__init__(FIAT.HsiehCloughTocher(cell))
+
+    def basis_transformation(self, coordinate_mapping):
+        # Jacobians at cell center
+        J = coordinate_mapping.jacobian_at([1/3, 1/3])
+
+        rns = coordinate_mapping.reference_normals()
+        pts = coordinate_mapping.physical_tangents()
+        pns = coordinate_mapping.physical_normals()
+
+        pel = coordinate_mapping.physical_edge_lengths()
+
+        d = self.cell.get_dimension()
+        ndof = self.space_dimension()
+        V = numpy.eye(ndof, dtype=object)
+        for multiindex in numpy.ndindex(V.shape):
+            V[multiindex] = Literal(V[multiindex])
+
+        voffset = d+1
+        for v in range(d+1):
+            s = voffset * v
+            for i in range(d):
+                for j in range(d):
+                    V[s+1+i, s+1+j] = J[j, i]
+
+        for e in range(3):
+            s = (d+1) * voffset + e
+            v0id, v1id = [i * voffset for i in range(3) if i != e]
+
+            nhat = partial_indexed(rns, (e, ))
+            t = partial_indexed(pts, (e, ))
+            n = partial_indexed(pns, (e, ))
+
+            Bn = J @ nhat / pel[e]
+            Bnn = Bn @ n
+            Bnt = Bn @ t
+
+            if self.avg:
+                Bnn = Bnn * pel[e]
+            V[s, s] = Bnn
+            V[s, v0id] = Literal(-1) * Bnt
+            V[s, v1id] = Bnt
+
+        # Patch up conditioning
+        h = coordinate_mapping.cell_size()
+        for v in range(d+1):
+            s = voffset * v
+            for k in range(d):
+                V[:, s+1+k] /= h[v]
+
+        for e in range(3):
+            v0id, v1id = [i for i in range(3) if i != e]
+            V[:, 9+e] *= 2 / (h[v0id] + h[v1id])
+        return ListTensor(V.T)
+
+
+class ReducedHsiehCloughTocher(PhysicallyMappedElement, ScalarFiatElement):
+    def __init__(self, cell, degree):
+        if degree != 3:
+            raise ValueError("Degree must be 3 for reduced HCT element")
+        if Citations is not None:
+            Citations().register("Clough1965")
+        super().__init__(FIAT.HsiehCloughTocher(cell, reduced=True))
+
+    def basis_transformation(self, coordinate_mapping):
+        # Jacobians at cell center
+        J = coordinate_mapping.jacobian_at([1/3, 1/3])
+
+        rns = coordinate_mapping.reference_normals()
+        pts = coordinate_mapping.physical_tangents()
+        # pns = coordinate_mapping.physical_normals()
+
+        pel = coordinate_mapping.physical_edge_lengths()
+
+        d = self.cell.get_dimension()
+        numbf = self._element.space_dimension()
+        ndof = self.space_dimension()
+        # rectangular to toss out the constraint dofs
+        V = numpy.eye(numbf, ndof, dtype=object)
+        for multiindex in numpy.ndindex(V.shape):
+            V[multiindex] = Literal(V[multiindex])
+
+        voffset = d+1
+        for v in range(d+1):
+            s = voffset * v
+            for i in range(d):
+                for j in range(d):
+                    V[s+1+i, s+1+j] = J[j, i]
+
+        for e in range(3):
+            s = (d+1) * voffset + e
+            v0id, v1id = [i * voffset for i in range(3) if i != e]
+
+            nhat = partial_indexed(rns, (e, ))
+            t = partial_indexed(pts, (e, ))
+
+            # n = partial_indexed(pns, (e, ))
+            # Bnn = (J @ nhat) @ n
+            # V[s, s] = Bnn
+
+            Bnt = (J @ nhat) @ t
+            V[s, v0id] = Literal(1/5) * Bnt / pel[e]
+            V[s, v1id] = Literal(-1) * V[s, v0id]
+
+            R = Literal(1/10) * Bnt * t
+            V[s, v0id + 1] = R[0]
+            V[s, v0id + 2] = R[1]
+            V[s, v1id + 1] = R[0]
+            V[s, v1id + 2] = R[1]
+
+        # Patch up conditioning
+        h = coordinate_mapping.cell_size()
+        for v in range(d+1):
+            s = voffset * v
+            for k in range(d):
+                V[:, s+1+k] /= h[v]
+        return ListTensor(V.T)
+
+    def entity_dofs(self):
+        edofs = deepcopy(super(ReducedHsiehCloughTocher, self).entity_dofs())
+        dim = 1
+        for entity in edofs[dim]:
+            edofs[dim][entity] = []
+        return edofs
+
+    @property
+    def index_shape(self):
+        return (9,)
+
+    def space_dimension(self):
+        return 9

finat/hdivcurl.py

@@ -26,6 +26,10 @@ def __init__(self, wrappee, transform):
     def cell(self):
         return self.wrappee.cell
 
+    @property
+    def complex(self):
+        return self.wrappee.complex
+
     @property
     def degree(self):
         return self.wrappee.degree

finat/mixed.py

@@ -37,6 +37,10 @@ def __init__(self, element, size, offset):
     def cell(self):
         return self.element.cell
 
+    @property
+    def complex(self):
+        return self.element.complex
+
     @property
     def degree(self):
         return self.element.degree

finat/physically_mapped.py

@@ -32,6 +32,15 @@
   eprint =       {1808.05513},
   primaryclass = {cs.MS}
 }""")
+    Citations().add("Clough1965", """
+@inproceedings{Clough1965,
+  author =       {R. W. Clough, J. L. Tocher},
+  title =        {Finite element stiffness matricess for analysis of plate bending},
+  booktitle =    {Proc. of the First Conf. on Matrix Methods in Struct. Mech},
+  year =         1965,
+  pages =        {515-546},
+}
+""")
     Citations().add("Argyris1968", """
 @Article{Argyris1968,
   author =       {J. H. Argyris and I. Fried and D. W. Scharpf},

finat/quadrature.py

@@ -1,16 +1,14 @@
 from abc import ABCMeta, abstractproperty
 from functools import reduce
 
-import numpy
-
 import gem
-from gem.utils import cached_property
-
-from FIAT.reference_element import LINE, QUADRILATERAL, TENSORPRODUCT
+import numpy
 from FIAT.quadrature import GaussLegendreQuadratureLineRule
 from FIAT.quadrature_schemes import create_quadrature as fiat_scheme
+from FIAT.reference_element import LINE, QUADRILATERAL, TENSORPRODUCT
+from gem.utils import cached_property
 
-from finat.point_set import PointSet, GaussLegendrePointSet, TensorPointSet
+from finat.point_set import GaussLegendrePointSet, PointSet, TensorPointSet
 
 
 def make_quadrature(ref_el, degree, scheme="default"):
@@ -44,7 +42,7 @@ def make_quadrature(ref_el, degree, scheme="default"):
     if degree < 0:
         raise ValueError("Need positive degree, not %d" % degree)
 
-    if ref_el.get_shape() == LINE:
+    if ref_el.get_shape() == LINE and not ref_el.is_macrocell():
         # FIAT uses Gauss-Legendre line quadature, however, since we
         # symbolically label it as such, we wish not to risk attaching
         # the wrong label in case FIAT changes.  So we explicitly ask