Merge pull request #11 from mscroggs/regge

Regge

Author: mscroggs

symfem/functionals.py

@@ -47,6 +47,37 @@ def entity_dim(self):
     name = "Point evaluation"
 
 
+class PointInnerProduct(BaseFunctional):
+    """An evaluation of an inner product at a point."""
+
+    def __init__(self, point, vec, entity_dim=None):
+        self.point = point
+        self.vec = vec
+        self._entity_dim = entity_dim
+
+    def eval(self, function):
+        """Apply to the functional to a function."""
+        v = subs(function, x, self.point)
+        tdim = len(self.vec)
+        return vdot(self.vec,
+                    tuple(vdot(v[tdim * i: tdim * (i + 1)], self.vec)
+                          for i in range(0, tdim)))
+
+    def dof_point(self):
+        """Get the location of the DOF in the cell."""
+        return self.point
+
+    def dof_direction(self):
+        """Get the location of the DOF in the cell."""
+        return self.vec
+
+    def entity_dim(self):
+        """Get the dimension of the entitiy this DOF is associated with."""
+        return self._entity_dim
+
+    name = "Point inner product"
+
+
 class DotPointEvaluation(BaseFunctional):
     """A point evaluation in a given direction."""
 

symfem/simplex.py

@@ -11,6 +11,7 @@
     TangentIntegralMoment,
     NormalIntegralMoment,
     IntegralMoment,
+    PointInnerProduct
 )
 
 
@@ -229,3 +230,38 @@ def __init__(self, reference, order):
 
     names = ["Crouzeix-Raviart", "CR"]
     min_order = 1
+
+
+class Regge(FiniteElement):
+    """A Regge element."""
+
+    def __init__(self, reference, order):
+        assert reference.name in ["triangle", "tetrahedron"]
+        if reference.tdim == 2:
+            poly = [(p[0], p[1], p[1], p[2])
+                    for p in polynomial_set(reference.tdim, 3, order)]
+        if reference.tdim == 3:
+            poly = [(p[0], p[1], p[3], p[1], p[2], p[4], p[3], p[4], p[5])
+                    for p in polynomial_set(reference.tdim, 6, order)]
+
+        dofs = []
+        for edim in range(1, 4):
+            for vs in reference.sub_entities(edim):
+                entity = reference.sub_entity_types[edim](
+                    vertices=tuple(reference.reference_vertices[i] for i in vs)
+                )
+                for i in product(range(1, order + 2), repeat=edim):
+                    if sum(i) < order + 2:
+                        for edge in entity.edges[::-1]:
+                            tangent = [b - a for a, b in zip(entity.vertices[edge[0]],
+                                                             entity.vertices[edge[1]])]
+                            dofs.append(PointInnerProduct(
+                                tuple(o + sum(sympy.Rational(a[j] * b, order + 2)
+                                              for a, b in zip(entity.axes, i[::-1]))
+                                      for j, o in enumerate(entity.origin)),
+                                tangent, entity_dim=edim))
+
+        super().__init__(reference, poly, dofs, reference.tdim, reference.tdim ** 2)
+
+    names = ["Regge"]
+    min_order = 0

test/test_against_basix.py

@@ -10,12 +10,14 @@
                  ("N2curl", "Nedelec 2nd kind H(curl)", range(1, 4)),
                  ("N1div", "Raviart-Thomas", range(1, 4)),
                  ("N2div", "Brezzi-Douglas-Marini", range(1, 4)),
+                 ("Regge", "Regge", range(0, 4)),
                  ("Crouzeix-Raviart", "Crouzeix-Raviart", [1])],
     "tetrahedron": [("P", "Lagrange", range(1, 4)), ("dP", "Discontinuous Lagrange", range(1, 4)),
                     ("N1curl", "Nedelec 1st kind H(curl)", range(1, 3)),
                     ("N2curl", "Nedelec 2nd kind H(curl)", range(1, 3)),
                     ("N1div", "Raviart-Thomas", range(1, 3)),
                     ("N2div", "Brezzi-Douglas-Marini", range(1, 3)),
+                    ("Regge", "Regge", range(0, 3)),
                     ("Crouzeix-Raviart", "Crouzeix-Raviart", [1])],
     "quadrilateral": [("Q", "Lagrange", range(1, 4)), ("dQ", "Discontinuous Lagrange", range(1, 4))],
     "hexahedron": [("Q", "Lagrange", range(1, 4)), ("dQ", "Discontinuous Lagrange", range(1, 4))]