Hermite on 1D manifolds

Author: pbrubeck

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,60 +9,41 @@
 
 
 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."""
+    """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):
-        entity_ids = {}
-        nodes = []
-        cur = 0
-
         # 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]):
+            cur = len(nodes)
             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))
+            if sd == 1:
+                # in 1D use normal derivative to support manifolds
+                nodes.append(functional.PointNormalDerivative(ref_el, v, verts[v]))
+            else:
+                nodes.extend(functional.PointDerivative(ref_el, verts[v], alpha)
+                             for alpha in polynomial_set.mis(sd, 1))
+            entity_ids[0][v].extend(range(cur, len(nodes)))
 
-            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] = []
-
+        # now only add dofs at the barycenter of the 2D entities
         if sd > 1:
             # face dof
             # point evaluation at barycenter
-            entity_ids[2] = {}
             for f in sorted(top[2]):
+                cur = len(nodes)
                 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] = []
+                nodes.append(functional.PointEvaluation(ref_el, pt))
+                entity_ids[2][f].extend(range(cur, len(nodes)))
 
         super().__init__(nodes, ref_el, entity_ids)
 

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,5 +1,6 @@
 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
@@ -15,18 +16,23 @@ 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", [
+                         finat.Hermite,
+                         ])
+def test_C1_interval(ref_to_phys, element):
+    check_zany_mapping(element, ref_to_phys[1.5])
+
+
 @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):