Merge pull request #21 from mscroggs/more-elements

More elements

Author: mscroggs

VERSION

@@ -1 +1 @@
-2021.1.1
+2021.1.2

setup.cfg

@@ -1,5 +1,5 @@
 [flake8]
-ignore = E501
+max-line-length = 100
 exclude = __init__.py
 min-version = 3.0
 

symfem/__init__.py

@@ -16,8 +16,8 @@
 with open(_os.path.join(_v_folder, "VERSION")) as _f:
     __version__ = _f.read().strip()
 
-_elementlist = {}
-
+_elementmap = {}
+_elementlist = []
 
 for _file in _os.listdir(_os.path.join(_folder, "elements")):
     if _file.endswith(".py") and "__init__" not in _file:
@@ -31,10 +31,11 @@
                 and issubclass(_element, _FiniteElement)
                 and _element != _FiniteElement
             ):
+                _elementlist.append(_element)
                 for _n in _element.names:
-                    if _n in _elementlist:
-                        assert _element == _elementlist[_n]
-                    _elementlist[_n] = _element
+                    if _n in _elementmap:
+                        assert _element == _elementmap[_n]
+                    _elementmap[_n] = _element
 
 
 def create_reference(cell_type, vertices=None):
@@ -89,13 +90,15 @@ def create_element(cell_type, element_type, order):
                           vector Lagrange, vP, Regge, Nedelec, Nedelec1, N1curl,
                           Nedelec2, N2curl, Raviart-Thomas, RT, N1div, dQ, NCE, RTCE,
                           Qcurl, Q, NCF, RTCF, Qdiv, vector Q, vQ,
-                          Brezzi-Douglas-Marini, BDM, N2div
+                          Brezzi-Douglas-Marini, BDM, N2div, Morley, Hermite,
+                          Mardal-Tai-Winther, MTW, Argyris
     order : int
         The order of the element.
     """
     reference = create_reference(cell_type)
 
-    if element_type in _elementlist:
-        return _elementlist[element_type](reference, order)
+    if element_type in _elementmap:
+        assert cell_type in _elementmap[element_type].references
+        return _elementmap[element_type](reference, order)
 
     raise ValueError(f"Unsupported element type: {element_type}")

symfem/core/calculus.py

@@ -0,0 +1,38 @@
+"""Functions to handle derivatives."""
+
+from .vectors import vdot
+from .symbolic import x, sym_sum
+
+
+def derivative(f, dir):
+    """Find the directional derivative of a function."""
+    return vdot(grad(f, len(dir)), dir)
+
+
+def grad(f, dim):
+    """Find the gradient of a scalar function."""
+    return tuple(f.diff(x[i]) for i in range(dim))
+
+
+def jacobian_component(f, component):
+    """Find a component of the Jacobian."""
+    return f.diff(x[component[0]]).diff(x[component[1]])
+
+
+def jacobian(f, dim):
+    """Find the Jacobian."""
+    return [[f.diff(x[i]).diff(x[j]) for i in range(dim)] for j in range(dim)]
+
+
+def div(f):
+    """Find the divergence of a vector function."""
+    return sym_sum(j.diff(x[i]) for i, j in enumerate(f))
+
+
+def curl(f):
+    """Find the curl of a 3D vector function."""
+    return (
+        f[2].diff(x[1]) - f[1].diff(x[2]),
+        f[0].diff(x[2]) - f[2].diff(x[0]),
+        f[1].diff(x[0]) - f[0].diff(x[1])
+    )

symfem/core/finite_element.py

@@ -118,14 +118,8 @@ def name(self):
 
 def make_integral_moment_dofs(
     reference,
-    vertices=None,
-    edges=None,
-    faces=None,
-    volumes=None,
-    cells=None,
-    facets=None,
-    ridges=None,
-    peaks=None,
+    vertices=None, edges=None, faces=None, volumes=None,
+    cells=None, facets=None, ridges=None, peaks=None
 ):
     """Generate DOFs due to integral moments on sub entities.
 

symfem/core/functionals.py

@@ -1,6 +1,7 @@
 """Functionals used to define the dual sets."""
 from .symbolic import subs, x, t
 from .vectors import vdot
+from .calculus import derivative, jacobian_component
 
 
 class BaseFunctional:
@@ -47,6 +48,66 @@ def entity_dim(self):
     name = "Point evaluation"
 
 
+class PointDirectionalDerivativeEvaluation(BaseFunctional):
+    """A point evaluation of a derivative in a fixed direction."""
+
+    def __init__(self, point, direction, entity_dim=None):
+        self.point = point
+        self.dir = direction
+        self._entity_dim = entity_dim
+
+    def eval(self, function):
+        """Apply to the functional to a function."""
+        return subs(derivative(function, self.dir), x, self.point)
+
+    def dof_point(self):
+        """Get the location of the DOF in the cell."""
+        return self.point
+
+    def dof_direction(self):
+        """Get the direction of the DOF."""
+        return self.dir
+
+    def entity_dim(self):
+        """Get the dimension of the entitiy this DOF is associated with."""
+        return self._entity_dim
+
+    name = "Point evaluation of directional derivative"
+
+
+class PointNormalDerivativeEvaluation(PointDirectionalDerivativeEvaluation):
+    """A point evaluation of a normal derivative."""
+
+    def __init__(self, point, edge):
+        super().__init__(point, edge.normal(), entity_dim=edge.tdim)
+        self.reference = edge
+
+    name = "Point evaluation of normal derivative"
+
+
+class PointComponentSecondDerivativeEvaluation(BaseFunctional):
+    """A point evaluation of a component of a second derivative."""
+
+    def __init__(self, point, component, entity_dim=None):
+        self.point = point
+        self.component = component
+        self._entity_dim = entity_dim
+
+    def eval(self, function):
+        """Apply to the functional to a function."""
+        return subs(jacobian_component(function, self.component), x, self.point)
+
+    def dof_point(self):
+        """Get the location of the DOF in the cell."""
+        return self.point
+
+    def entity_dim(self):
+        """Get the dimension of the entitiy this DOF is associated with."""
+        return self._entity_dim
+
+    name = "Point evaluation of Jacobian component"
+
+
 class PointInnerProduct(BaseFunctional):
     """An evaluation of an inner product at a point."""
 

symfem/core/symbolic.py

@@ -22,3 +22,11 @@ def subs(f, vars, values):
         return (
             f.subs(vars[0], values[0]).subs(vars[1], values[1]).subs(vars[2], values[2])
         )
+
+
+def sym_sum(ls):
+    """Symbolically computes the sum of a list."""
+    out = zero
+    for i in ls:
+        out += i
+    return out

symfem/core/vectors.py

@@ -11,6 +11,14 @@ def vsub(v, w):
         return v - w
 
 
+def vadd(v, w):
+    """Add two vectors."""
+    try:
+        return tuple(i + j for i, j in zip(v, w))
+    except TypeError:
+        return v + w
+
+
 def vdiv(v, a):
     """Divide a vector by a scalar."""
     try:

symfem/elements/argyris.py

@@ -0,0 +1,43 @@
+"""Argyris elements on simplices."""
+
+from ..core.finite_element import FiniteElement
+from ..core.polynomials import polynomial_set
+from ..core.functionals import (PointEvaluation, PointDirectionalDerivativeEvaluation,
+                                PointNormalDerivativeEvaluation,
+                                PointComponentSecondDerivativeEvaluation)
+from ..core.symbolic import sym_sum
+
+
+class Argyris(FiniteElement):
+    """Argyris finite element."""
+
+    def __init__(self, reference, order):
+        from symfem import create_reference
+        assert order == 5
+        assert reference.name == "triangle"
+        dofs = []
+        for vs in reference.sub_entities(0):
+            dofs.append(PointEvaluation(vs, entity_dim=0))
+            for i in range(reference.tdim):
+                dir = tuple(1 if i == j else 0 for j in range(reference.tdim))
+                dofs.append(PointDirectionalDerivativeEvaluation(vs, dir, entity_dim=0))
+            for i in range(reference.tdim):
+                for j in range(i + 1):
+                    dofs.append(PointComponentSecondDerivativeEvaluation(
+                        vs, (i, j), entity_dim=0))
+        for vs in reference.sub_entities(1):
+            sub_ref = create_reference(
+                reference.sub_entity_types[1],
+                vertices=[reference.reference_vertices[v] for v in vs])
+            midpoint = tuple(sym_sum(i) / len(i)
+                             for i in zip(*[reference.vertices[i] for i in vs]))
+            dofs.append(PointNormalDerivativeEvaluation(midpoint, sub_ref))
+
+        super().__init__(
+            reference, polynomial_set(reference.tdim, 1, order), dofs, reference.tdim, 1
+        )
+
+    names = ["Argyris"]
+    references = ["triangle"]
+    min_order = 5
+    max_order = 5

symfem/elements/bdm.py

@@ -22,4 +22,5 @@ def __init__(self, reference, order):
         super().__init__(reference, poly, dofs, reference.tdim, reference.tdim)
 
     names = ["Brezzi-Douglas-Marini", "BDM", "N2div"]
+    references = ["triangle", "tetrahedron"]
     min_order = 1