Add Direct serendipity element (#83)

* add CiarletElement class

* pydocstyle

* Add Direct serendipity element

* flake

Author: mscroggs

demo/custom_element.py

@@ -1,11 +1,11 @@
 import symfem
 import sympy
-from symfem.core.finite_element import FiniteElement
+from symfem.core.finite_element import CiarletElement
 from symfem.core.symbolic import x
 from symfem.core.functionals import PointEvaluation
 
 
-class CustomElement(FiniteElement):
+class CustomElement(CiarletElement):
     """Custom element on a quadrilateral."""
 
     def __init__(self, reference, order, variant):

symfem/core/finite_element.py

@@ -9,74 +9,22 @@
 class FiniteElement:
     """Abstract finite element."""
 
-    def __init__(self, reference, order, basis, dofs, domain_dim, range_dim,
+    def __init__(self, reference, order, space_dim, domain_dim, range_dim,
                  range_shape=None):
-        assert len(basis) == len(dofs)
         self.reference = reference
         self.order = order
-        self.basis = basis
-        self.dofs = dofs
+        self.space_dim = space_dim
         self.domain_dim = domain_dim
         self.range_dim = range_dim
         self.range_shape = range_shape
-        self.space_dim = len(dofs)
-        self._basis_functions = None
-        self._reshaped_basis_functions = None
 
     def entity_dofs(self, entity_dim, entity_number):
         """Get the numbers of the DOFs associated with the given entity."""
-        return [i for i, j in enumerate(self.dofs) if j.entity == (entity_dim, entity_number)]
-
-    def get_polynomial_basis(self, reshape=True):
-        """Get the polynomial basis for the element."""
-        if reshape and self.range_shape is not None:
-            if len(self.range_shape) != 2:
-                raise NotImplementedError
-            assert self.range_shape[0] * self.range_shape[1] == self.range_dim
-            return [sympy.Matrix(
-                [b[i * self.range_shape[1]: (i + 1) * self.range_shape[1]]
-                 for i in range(self.range_shape[0])]) for b in self.basis]
-
-        return self.basis
-
-    def get_dual_matrix(self):
-        """Get the dual matrix."""
-        mat = []
-        for b in self.basis:
-            row = []
-            for d in self.dofs:
-                row.append(d.eval(b))
-            mat.append(row)
-        return sympy.Matrix(mat)
+        raise NotImplementedError()
 
     def get_basis_functions(self, reshape=True):
         """Get the basis functions of the element."""
-        if self._basis_functions is None:
-            minv = self.get_dual_matrix().inv("LU")
-            self._basis_functions = []
-            if self.range_dim == 1:
-                # Scalar space
-                for i, dof in enumerate(self.dofs):
-                    self._basis_functions.append(
-                        sym_sum(c * d for c, d in zip(minv.row(i), self.basis)))
-            else:
-                # Vector space
-                for i, dof in enumerate(self.dofs):
-                    b = [zero for i in self.basis[0]]
-                    for c, d in zip(minv.row(i), self.basis):
-                        for j, d_j in enumerate(d):
-                            b[j] += c * d_j
-                    self._basis_functions.append(b)
-
-        if reshape and self.range_shape is not None:
-            if len(self.range_shape) != 2:
-                raise NotImplementedError
-            assert self.range_shape[0] * self.range_shape[1] == self.range_dim
-            return [sympy.Matrix(
-                [b[i * self.range_shape[1]: (i + 1) * self.range_shape[1]]
-                 for i in range(self.range_shape[0])]) for b in self._basis_functions]
-
-        return self._basis_functions
+        raise NotImplementedError()
 
     def get_basis_function(self, n):
         """Get a single basis function of the element."""
@@ -147,3 +95,103 @@ def name(self):
         return self.names[0]
 
     names = []
+
+
+class CiarletElement(FiniteElement):
+    """Finite element defined using the Ciarlet definition."""
+
+    def __init__(self, reference, order, basis, dofs, domain_dim, range_dim,
+                 range_shape=None):
+        super().__init__(reference, order, len(dofs), domain_dim, range_dim, range_shape)
+        assert len(basis) == len(dofs)
+        self.basis = basis
+        self.dofs = dofs
+        self._basis_functions = None
+        self._reshaped_basis_functions = None
+
+    def entity_dofs(self, entity_dim, entity_number):
+        """Get the numbers of the DOFs associated with the given entity."""
+        return [i for i, j in enumerate(self.dofs) if j.entity == (entity_dim, entity_number)]
+
+    def get_polynomial_basis(self, reshape=True):
+        """Get the polynomial basis for the element."""
+        if reshape and self.range_shape is not None:
+            if len(self.range_shape) != 2:
+                raise NotImplementedError
+            assert self.range_shape[0] * self.range_shape[1] == self.range_dim
+            return [sympy.Matrix(
+                [b[i * self.range_shape[1]: (i + 1) * self.range_shape[1]]
+                 for i in range(self.range_shape[0])]) for b in self.basis]
+
+        return self.basis
+
+    def get_dual_matrix(self):
+        """Get the dual matrix."""
+        mat = []
+        for b in self.basis:
+            row = []
+            for d in self.dofs:
+                row.append(d.eval(b))
+            mat.append(row)
+        return sympy.Matrix(mat)
+
+    def get_basis_functions(self, reshape=True):
+        """Get the basis functions of the element."""
+        if self._basis_functions is None:
+            minv = self.get_dual_matrix().inv("LU")
+            self._basis_functions = []
+            if self.range_dim == 1:
+                # Scalar space
+                for i, dof in enumerate(self.dofs):
+                    self._basis_functions.append(
+                        sym_sum(c * d for c, d in zip(minv.row(i), self.basis)))
+            else:
+                # Vector space
+                for i, dof in enumerate(self.dofs):
+                    b = [zero for i in self.basis[0]]
+                    for c, d in zip(minv.row(i), self.basis):
+                        for j, d_j in enumerate(d):
+                            b[j] += c * d_j
+                    self._basis_functions.append(b)
+
+        if reshape and self.range_shape is not None:
+            if len(self.range_shape) != 2:
+                raise NotImplementedError
+            assert self.range_shape[0] * self.range_shape[1] == self.range_dim
+            return [sympy.Matrix(
+                [b[i * self.range_shape[1]: (i + 1) * self.range_shape[1]]
+                 for i in range(self.range_shape[0])]) for b in self._basis_functions]
+
+        return self._basis_functions
+
+    names = []
+
+
+class DirectElement(FiniteElement):
+    """Finite element defined directly."""
+
+    def __init__(self, reference, order, basis_functions, basis_entities, domain_dim, range_dim,
+                 range_shape=None):
+        super().__init__(reference, order, len(basis_functions), domain_dim, range_dim,
+                         range_shape)
+        self._basis_entities = basis_entities
+        self._basis_functions = basis_functions
+        self._reshaped_basis_functions = None
+
+    def entity_dofs(self, entity_dim, entity_number):
+        """Get the numbers of the DOFs associated with the given entity."""
+        return [i for i, j in enumerate(self._basis_entities) if j == (entity_dim, entity_number)]
+
+    def get_basis_functions(self, reshape=True):
+        """Get the basis functions of the element."""
+        if reshape and self.range_shape is not None:
+            if len(self.range_shape) != 2:
+                raise NotImplementedError
+            assert self.range_shape[0] * self.range_shape[1] == self.range_dim
+            return [sympy.Matrix(
+                [b[i * self.range_shape[1]: (i + 1) * self.range_shape[1]]
+                 for i in range(self.range_shape[0])]) for b in self._basis_functions]
+
+        return self._basis_functions
+
+    names = []

symfem/create.py

@@ -142,6 +142,7 @@ def create_element(cell_type, element_type, order, variant="equispaced"):
         serendipity, S,
         serendipity Hcurl, Scurl, BDMCE, AAE,
         serendipity Hdiv, Sdiv, BDMCF, AAF,
+        direct serendipity,
         Regge,
         Nedelec, Nedelec1, N1curl,
         Nedelec2, N2curl,

symfem/elements/argyris.py

@@ -4,15 +4,15 @@
 (Arygris, Fried, Scharpf, 1968)
 """
 
-from ..core.finite_element import FiniteElement
+from ..core.finite_element import CiarletElement
 from ..core.polynomials import polynomial_set
 from ..core.functionals import (PointEvaluation, PointDirectionalDerivativeEvaluation,
                                 PointNormalDerivativeEvaluation,
                                 PointComponentSecondDerivativeEvaluation)
 from ..core.symbolic import sym_sum
 
 
-class Argyris(FiniteElement):
+class Argyris(CiarletElement):
     """Argyris finite element."""
 
     def __init__(self, reference, order, variant):

symfem/elements/aw.py

@@ -5,15 +5,15 @@
 """
 
 import sympy
-from ..core.finite_element import FiniteElement
+from ..core.finite_element import CiarletElement
 from ..core.polynomials import polynomial_set
 from ..core.functionals import (PointInnerProduct, InnerProductIntegralMoment,
                                 VecIntegralMoment, IntegralMoment)
 from ..core.symbolic import zero, one, x
 from .lagrange import DiscontinuousLagrange
 
 
-class ArnoldWinther(FiniteElement):
+class ArnoldWinther(CiarletElement):
     """An Arnold-Winther element."""
 
     def __init__(self, reference, order, variant):

symfem/elements/bddm.py

@@ -4,7 +4,7 @@
 (Brezzi, Douglas, Duran, Fortin, 1987)
 """
 
-from ..core.finite_element import FiniteElement
+from ..core.finite_element import CiarletElement
 from ..core.moments import make_integral_moment_dofs
 from ..core.polynomials import polynomial_set
 from ..core.symbolic import x, zero
@@ -30,7 +30,7 @@ def bddf_polyset(reference, order):
     return pset
 
 
-class BDDF(FiniteElement):
+class BDDF(CiarletElement):
     """Brezzi-Douglas-Duran-Fortin Hdiv finite element."""
 
     def __init__(self, reference, order, variant):

symfem/elements/bdfm.py

@@ -4,7 +4,7 @@
 (Brezzi, Douglas, Fortin, Marini, 1987)
 """
 
-from ..core.finite_element import FiniteElement
+from ..core.finite_element import CiarletElement
 from ..core.moments import make_integral_moment_dofs
 from ..core.polynomials import polynomial_set
 from ..core.symbolic import x, zero
@@ -44,7 +44,7 @@ def bdfm_polyset(reference, order):
     return pset
 
 
-class BDFM(FiniteElement):
+class BDFM(CiarletElement):
     """Brezzi-Douglas-Fortin-Marini Hdiv finite element."""
 
     def __init__(self, reference, order, variant):

symfem/elements/bdm.py

@@ -4,15 +4,15 @@
 (Brezzi, Douglas, Marini, 1985)
 """
 
-from ..core.finite_element import FiniteElement
+from ..core.finite_element import CiarletElement
 from ..core.moments import make_integral_moment_dofs
 from ..core.polynomials import polynomial_set
 from ..core.functionals import NormalIntegralMoment, IntegralMoment
 from .lagrange import DiscontinuousLagrange
 from .nedelec import NedelecFirstKind
 
 
-class BDM(FiniteElement):
+class BDM(CiarletElement):
     """Brezzi-Douglas-Marini Hdiv finite element."""
 
     def __init__(self, reference, order, variant):

symfem/elements/bell.py

@@ -3,15 +3,15 @@
 This element's definition is given in https://doi.org/10.1002/nme.1620010108 (Bell, 1969)
 """
 
-from ..core.finite_element import FiniteElement
+from ..core.finite_element import CiarletElement
 from ..core.moments import make_integral_moment_dofs
 from ..core.polynomials import polynomial_set
 from ..core.functionals import (PointEvaluation, NormalDerivativeIntegralMoment,
                                 DerivativePointEvaluation)
 from .lagrange import DiscontinuousLagrange
 
 
-class Bell(FiniteElement):
+class Bell(CiarletElement):
     """Bell finite element."""
 
     def __init__(self, reference, order, variant):

symfem/elements/bernardi_raugel.py

@@ -5,14 +5,14 @@
 """
 
 from ..core import mappings
-from ..core.finite_element import FiniteElement
+from ..core.finite_element import CiarletElement
 from ..core.moments import make_integral_moment_dofs
 from ..core.polynomials import polynomial_set
 from ..core.functionals import NormalIntegralMoment, DotPointEvaluation
 from .lagrange import Lagrange, DiscontinuousLagrange
 
 
-class BernardiRaugel(FiniteElement):
+class BernardiRaugel(CiarletElement):
     """Bernardi-Raugel Hdiv finite element."""
 
     def __init__(self, reference, order, variant):