HCT: basis transformation for hierarchical edge functions (#131)

Author: pbrubeck

finat/argyris.py

@@ -9,18 +9,28 @@
 from finat.physically_mapped import PhysicallyMappedElement, Citations
 
 
-def _edge_transform(V, voffset, fiat_cell, moment_deg, coordinate_mapping, avg=False):
-    """Basis transformation for integral edge moments."""
-
-    J = coordinate_mapping.jacobian_at([1/3, 1/3])
+def _edge_transform(V, vorder, eorder, fiat_cell, coordinate_mapping, avg=False):
+    """Basis transformation for integral edge moments.
+
+    :arg V: the transpose of the basis transformation.
+    :arg vorder: the jet order at vertices, matching the Jacobi weights in the
+                 normal derivative moments on edges.
+    :arg eorder: the order of the normal derivative moments.
+    :arg fiat_cell: the reference triangle.
+    :arg coordinate_mapping: the coordinate mapping.
+    :kwarg avg: are we scaling integrals by dividing by the edge length?
+    """
+    sd = fiat_cell.get_spatial_dimension()
+    J = coordinate_mapping.jacobian_at(fiat_cell.make_points(sd, 0, sd+1)[0])
     rns = coordinate_mapping.reference_normals()
     pts = coordinate_mapping.physical_tangents()
     pns = coordinate_mapping.physical_normals()
     pel = coordinate_mapping.physical_edge_lengths()
 
+    # number of DOFs per vertex/edge
+    voffset = comb(sd + vorder, vorder)
+    eoffset = 2 * eorder + 1
     top = fiat_cell.get_topology()
-    eoffset = 2 * moment_deg - 1
-    toffset = moment_deg - 1
     for e in sorted(top[1]):
         nhat = partial_indexed(rns, (e, ))
         n = partial_indexed(pns, (e, ))
@@ -33,15 +43,16 @@ def _edge_transform(V, voffset, fiat_cell, moment_deg, coordinate_mapping, avg=F
 
         v0id, v1id = (v * voffset for v in top[1][e])
         s0 = len(top[0]) * voffset + e * eoffset
-        for k in range(moment_deg):
+        for k in range(eorder+1):
             s = s0 + k
-            P1 = Literal(comb(k + 2, k))
+            # Jacobi polynomial at the endpoints
+            P1 = comb(k + vorder, k)
             P0 = -(-1)**k * P1
             V[s, s] = Bnn
             V[s, v1id] = P1 * Bnt
             V[s, v0id] = P0 * Bnt
             if k > 0:
-                V[s, s + toffset] = -1 * Bnt
+                V[s, s + eorder] = -1 * Bnt
 
 
 class Argyris(PhysicallyMappedElement, ScalarFiatElement):
@@ -68,7 +79,8 @@ def basis_transformation(self, coordinate_mapping):
 
         sd = self.cell.get_spatial_dimension()
         top = self.cell.get_topology()
-        voffset = (sd+1)*sd//2 + sd + 1
+        vorder = 2
+        voffset = comb(sd + vorder, vorder)
         for v in sorted(top[0]):
             s = voffset * v
             for i in range(sd):
@@ -84,9 +96,9 @@ def basis_transformation(self, coordinate_mapping):
             V[s+5, s+4] = 2*J[0, 1]*J[1, 1]
             V[s+5, s+5] = J[1, 1]*J[1, 1]
 
-        q = self.degree - 4
+        eorder = self.degree - 5
         if self.variant == "integral":
-            _edge_transform(V, voffset, self.cell, q, coordinate_mapping, avg=self.avg)
+            _edge_transform(V, vorder, eorder, self.cell, coordinate_mapping, avg=self.avg)
         else:
             rns = coordinate_mapping.reference_normals()
             pns = coordinate_mapping.physical_normals()
@@ -128,10 +140,10 @@ def basis_transformation(self, coordinate_mapping):
                 V[:, voffset*v+3+k] *= 1 / (h[v]*h[v])
 
         if self.variant == "point":
-            eoffset = 2 * q - 1
+            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+q] *= 2 / (h[v0] + h[v1])
+                V[:, s:s+eorder+1] *= 2 / (h[v0] + h[v1])
 
         return ListTensor(V.T)

finat/hct.py

@@ -30,15 +30,16 @@ def basis_transformation(self, coordinate_mapping):
 
         sd = self.cell.get_dimension()
         top = self.cell.get_topology()
-        voffset = sd + 1
+        voffset = 1 + sd
         for v in sorted(top[0]):
             s = voffset * v
             for i in range(sd):
                 for j in range(sd):
                     V[s+1+i, s+1+j] = J[j, i]
 
-        q = self.degree - 2
-        _edge_transform(V, voffset, self.cell, q, coordinate_mapping, avg=self.avg)
+        vorder = 1
+        eorder = self.degree - 3
+        _edge_transform(V, vorder, eorder, self.cell, coordinate_mapping, avg=self.avg)
 
         # Patch up conditioning
         h = coordinate_mapping.cell_size()