Higher-order Hermite on 1D manifolds

FIAT/expansions.py

@@ -580,7 +580,7 @@ def __init__(self, ref_el, **kwargs):
     def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         """Returns a dict of tabulations such that
         tabulations[alpha][i, j] = D^alpha phi_i(pts[j])."""
-        if self.variant is not None:
+        if self.variant is not None or self.ref_el.get_spatial_dimension() != 1:
             return super()._tabulate_on_cell(n, pts, order=order, cell=cell, direction=direction)
 
         A, b = self.affine_mappings[cell]

FIAT/hermite.py

@@ -9,70 +9,58 @@
 
 
 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
+    """The dual basis for Hermite elements.  This class works for
+    simplices of any dimension.  Nodes are first order jet at
+    vertices and point evaluation at barycenters of 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()
-
-        # get jet at each vertex
+        sd = ref_el.get_topological_dimension()
+        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
+        nodes = []
 
-        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))
+            if sd == 1:
+                # use normal derivative to support manifolds in 1D
+                nodes.append(functional.PointNormalDerivative(ref_el, v, pt))
+            else:
+                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:
+            assert degree == 3
+            # 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."""
 
-    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 variant is None:
+            variant = "gll"
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
+        dual = CubicHermiteDualSet(ref_el, degree, variant=variant)
 
-        super().__init__(poly_set, dual, 3)
+        super().__init__(poly_set, dual, degree)

FIAT/reference_element.py

@@ -1001,6 +1001,11 @@ def __init__(self):
                     1: edges}
         super().__init__(LINE, verts, topology)
 
+    def compute_normal(self, i):
+        "UFC consistent normal"
+        n = self.compute_tangents(1, 0)[0]
+        return n / numpy.linalg.norm(n)
+
 
 class DefaultTriangle(DefaultSimplex):
     """This is the reference triangle with vertices (-1.0,-1.0),

finat/hermite.py

@@ -1,32 +1,38 @@
 import FIAT
-from gem import ListTensor
+import numpy
+from gem import ListTensor, partial_indexed
 
 from finat.citations import cite
 from finat.fiat_elements import ScalarFiatElement
 from finat.physically_mapped import identity, PhysicallyMappedElement
 
 
 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.CubicHermite(cell, degree=degree, variant=variant))
 
     def basis_transformation(self, coordinate_mapping):
         Js = [coordinate_mapping.jacobian_at(vertex)
               for vertex in self.cell.get_vertices()]
 
+        pns = coordinate_mapping.physical_normals()
         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
+
+            gdim, tdim = J.shape
+            if gdim != tdim:
+                assert tdim == 1
+                J = partial_indexed(pns, (i,)) @ J
+
+            Jnp = numpy.reshape([J[i] for i in numpy.ndindex(J.shape)], J.shape)
+            M[numpy.ix_(vids, vids)] = Jnp * (1 / h[i])
 
         return ListTensor(M)

test/finat/conftest.py

@@ -107,18 +107,19 @@ 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(int(dim)) for dim in (2, 2.5, 3)}
+    K = {dim: FIAT.ufc_simplex(int(dim)) for dim in (1.5, 2, 2.5, 3)}
     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),
                      (0.01, 0.987, -.23),
                      (-0.1, -0.2, 1.38))
+    K[1.5].vertices = K[2].vertices[:-1]
     K[2.5].vertices = K[3].vertices[:-1]
     return K
 

test/finat/test_zany_mapping.py

@@ -119,6 +119,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.5], degree)
+
+
 @pytest.mark.parametrize("element", [
                          finat.Morley,
                          finat.Hermite,
@@ -132,6 +139,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):