Higher-order Hermite on 1D manifolds

FIAT/__init__.py

@@ -39,7 +39,7 @@
 from FIAT.serendipity import Serendipity
 from FIAT.brezzi_douglas_marini_cube import BrezziDouglasMariniCubeEdge, BrezziDouglasMariniCubeFace
 from FIAT.discontinuous_pc import DPC
-from FIAT.hermite import CubicHermite
+from FIAT.hermite import Hermite
 from FIAT.lagrange import Lagrange
 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre
 from FIAT.gauss_legendre import GaussLegendre
@@ -89,7 +89,7 @@
                       "DPC": DPC,
                       "Discontinuous Taylor": DiscontinuousTaylor,
                       "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas,
-                      "Hermite": CubicHermite,
+                      "Hermite": Hermite,
                       "Nonconforming Wu-Xu": WuXuH3NC,
                       "Nonconforming Robust Wu-Xu": WuXuRobustH3NC,
                       "Hsieh-Clough-Tocher": HsiehCloughTocher,

FIAT/hermite.py

@@ -8,71 +8,59 @@
 from FIAT import finite_element, polynomial_set, dual_set, functional
 
 
-class CubicHermiteDualSet(dual_set.DualSet):
-    """The dual basis for Lagrange elements.  This class works for
-    simplices of any dimension.  Nodes are point evaluation at
-    equispaced points."""
-
-    def __init__(self, ref_el):
-        entity_ids = {}
-        nodes = []
-        cur = 0
-
+class HermiteDualSet(dual_set.DualSet):
+    """The Hermite dual set is defined in 1D for any degree, with degrees of
+    freedom given by the first-order jet at the vertices and point
+    evaluations at interior points of the interval.  In higher dimensions,
+    it is defined only for degree 3, with the first-order jet at the vertices
+    and point evaluations at the barycenters of the 2D entities."""
+
+    def __init__(self, ref_el, degree, variant=None):
         # make nodes by getting points
         # need to do this dimension-by-dimension, facet-by-facet
         top = ref_el.get_topology()
-        verts = ref_el.get_vertices()
         sd = ref_el.get_spatial_dimension()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
-        # get jet at each vertex
-
-        entity_ids[0] = {}
+        # get first order jet at each vertex
         for v in sorted(top[0]):
-            nodes.append(functional.PointEvaluation(ref_el, verts[v]))
-            pd = functional.PointDerivative
-            for i in range(sd):
-                alpha = [0] * sd
-                alpha[i] = 1
-
-                nodes.append(pd(ref_el, verts[v], alpha))
-
-            entity_ids[0][v] = list(range(cur, cur + 1 + sd))
-            cur += sd + 1
-
-        # now only have dofs at the barycenter, which is the
-        # maximal dimension
-        # no edge dof
-
-        entity_ids[1] = {}
-        for i in top[1]:
-            entity_ids
-            entity_ids[1][i] = []
-
-        if sd > 1:
-            # face dof
-            # point evaluation at barycenter
-            entity_ids[2] = {}
+            pt, = ref_el.make_points(0, v, degree, variant=variant)
+            cur = len(nodes)
+            nodes.append(functional.PointEvaluation(ref_el, pt))
+            nodes.extend(functional.PointDerivative(ref_el, pt, alpha)
+                         for alpha in polynomial_set.mis(sd, 1))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
+
+        if sd == 1:
+            # edge dofs: point evaluations to support higher order in 1D
+            for e in sorted(top[1]):
+                cur = len(nodes)
+                pts = ref_el.make_points(1, e, degree-2, variant=variant)
+                nodes.extend(functional.PointEvaluation(ref_el, pt) for pt in pts)
+                entity_ids[1][e].extend(range(cur, len(nodes)))
+        else:
+            # no edge dof
+            # face dof: point evaluation at barycenter
             for f in sorted(top[2]):
-                pt = ref_el.make_points(2, f, 3)[0]
-                n = functional.PointEvaluation(ref_el, pt)
-                nodes.append(n)
-                entity_ids[2][f] = list(range(cur, cur + 1))
-                cur += 1
-
-            for dim in range(3, sd + 1):
-                entity_ids[dim] = {}
-                for facet in top[dim]:
-                    entity_ids[dim][facet] = []
+                cur = len(nodes)
+                pt, = ref_el.make_points(2, f, degree, variant=variant)
+                nodes.append(functional.PointEvaluation(ref_el, pt))
+                entity_ids[2][f].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 
 
-class CubicHermite(finite_element.CiarletElement):
-    """The cubic Hermite finite element.  It is what it is."""
+class Hermite(finite_element.CiarletElement):
+    """The Hermite finite element.  It has any degree of at least three
+    on intervals, and degree three on higher dimensional simplices."""
 
-    def __init__(self, ref_el, deg=3):
-        assert deg == 3
-        poly_set = polynomial_set.ONPolynomialSet(ref_el, 3)
-        dual = CubicHermiteDualSet(ref_el)
+    def __init__(self, ref_el, degree=3, variant=None):
+        if ref_el.get_spatial_dimension() > 1 and degree != 3:
+            raise ValueError("Hermite elements in more than one dimension must have degree 3.")
+        if variant is None:
+            variant = "gll"
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
+        dual = HermiteDualSet(ref_el, degree, variant=variant)
 
-        super().__init__(poly_set, dual, 3)
+        super().__init__(poly_set, dual, degree)

finat/hermite.py

@@ -1,4 +1,5 @@
 import FIAT
+import numpy
 from gem import ListTensor
 
 from finat.citations import cite
@@ -7,26 +8,29 @@
 
 
 class Hermite(PhysicallyMappedElement, ScalarFiatElement):
-    def __init__(self, cell, degree=3):
+    def __init__(self, cell, degree=3, variant=None):
         cite("Ciarlet1972")
-        super().__init__(FIAT.CubicHermite(cell))
+        super().__init__(FIAT.Hermite(cell, degree=degree, variant=variant))
 
     def basis_transformation(self, coordinate_mapping):
-        Js = [coordinate_mapping.jacobian_at(vertex)
-              for vertex in self.cell.get_vertices()]
+        vertices = self.cell.get_vertices()
+        if self.cell.get_spatial_dimension() == 1:
+            # The derivative along the tangent maps by the signed Jacobian
+            # determinant, which carries the cell orientation on manifolds.
+            Js = [ListTensor([[coordinate_mapping.detJ_at(vertex)]]) for vertex in vertices]
+        else:
+            Js = [coordinate_mapping.jacobian_at(vertex) for vertex in vertices]
 
         h = coordinate_mapping.cell_size()
 
-        d = self.cell.get_dimension()
         M = identity(self.space_dimension())
 
-        cur = 0
-        for i in range(d+1):
-            cur += 1  # skip the vertex
+        entity_ids = self.entity_dofs()
+        for i in entity_ids[0]:
+            # skip the PointEvaluation DOF
+            vids = entity_ids[0][i][1:]
             J = Js[i]
-            for j in range(d):
-                for k in range(d):
-                    M[cur+j, cur+k] = J[j, k] / h[i]
-            cur += d
+            Jnp = numpy.reshape([J[k] for k in numpy.ndindex(J.shape)], J.shape)
+            M[numpy.ix_(vids, vids)] = Jnp * (1 / h[i])
 
         return ListTensor(M)

test/FIAT/unit/test_fiat.py

@@ -50,7 +50,7 @@
 from FIAT.hu_zhang import HuZhang                               # noqa: F401
 from FIAT.bernardi_raugel import BernardiRaugel                 # noqa: F401
 from FIAT.argyris import Argyris                                # noqa: F401
-from FIAT.hermite import CubicHermite                           # noqa: F401
+from FIAT.hermite import Hermite                                # noqa: F401
 from FIAT.morley import Morley                                  # noqa: F401
 from FIAT.hct import HsiehCloughTocher                          # noqa: F401
 from FIAT.c2_elements import AlfeldC2, BrambleZlamalC2          # noqa: F401
@@ -353,9 +353,10 @@ def __init__(self, a, b):
     "Argyris(T, 6, 'integral')",
     "WuXuH3NC(T, 4)",
     "WuXuRobustH3NC(T, 7)",
-    "CubicHermite(I)",
-    "CubicHermite(T)",
-    "CubicHermite(S)",
+    "Hermite(I)",
+    "Hermite(I, 4)",
+    "Hermite(T)",
+    "Hermite(S)",
     "Morley(T)",
     "Morley(S)",
     "BernardiRaugel(T)",

test/FIAT/unit/test_pointwise_dual.py

@@ -8,7 +8,7 @@
 import numpy
 
 from FIAT import (
-    BrezziDouglasMarini, Morley, Argyris, CubicHermite)
+    BrezziDouglasMarini, Morley, Argyris, Hermite)
 
 from FIAT.reference_element import (
     UFCTriangle,
@@ -20,7 +20,7 @@
 
 
 @pytest.mark.parametrize("element",
-                         [CubicHermite(T),
+                         [Hermite(T),
                           Morley(T),
                           Argyris(T),
                           BrezziDouglasMarini(T, 1, variant="integral")])

test/finat/conftest.py

@@ -102,13 +102,14 @@ def scaled_simplex(dim, scale):
 
 @pytest.fixture
 def ref_el():
-    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
+    K = {dim: FIAT.ufc_simplex(dim) for dim in (1, 2, 3)}
     return K
 
 
 @pytest.fixture
 def phys_el():
-    K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)}
+    K = {dim: FIAT.ufc_simplex(dim) for dim in (1, 2, 3)}
+    K[1].vertices = ((0.1,), (1.27,))
     K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84))
     K[3].vertices = ((0, 0, 0),
                      (1., 0.1, -0.37),

test/finat/test_zany_mapping.py

@@ -103,6 +103,13 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs):
     assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), pp.pformat((np.round(error, 8).tolist(), *inds))
 
 
+@pytest.mark.parametrize("element, degree", [
+    *((finat.Hermite, k) for k in range(3, 6)),
+])
+def test_C1_interval(ref_to_phys, element, degree):
+    check_zany_mapping(element, ref_to_phys[1], degree)
+
+
 @pytest.mark.parametrize("element", [
                          finat.Morley,
                          finat.Hermite,
@@ -116,6 +123,7 @@ def test_C1_triangle(ref_to_phys, element):
 
 @pytest.mark.parametrize("element", [
                          finat.Morley,
+                         finat.Hermite,
                          finat.Walkington,
                          ])
 def test_C1_tetrahedron(ref_to_phys, element):