Legendre and Lobatto elements (#231)

* Implement Legendre and Lobatto Lagrange elements

* isort

* flake8

* pydocstyle

* mypy

* flake

* isort

* setup.py

* fix tests

* remove print

* fix DPc element

* cache date?

* isort

* fix tests

Author: mscroggs

CHANGELOG_SINCE_LAST_VERSION.md

@@ -5,3 +5,4 @@
 - Corrected C1 and higher order C tests
 - Allow element creation on non-default references
 - Corrected Bell element
+- Corrected legendre and lobatto variants of Lagrange

setup.py

@@ -34,7 +34,7 @@
         author_email="[email protected]",
         maintainer_email="[email protected]",
         url="https://github.com/mscroggs/symfem",
-        packages=["symfem", "symfem.elements"],
+        packages=["symfem", "symfem.elements", "symfem.polynomials"],
         package_data={"symfem": ["py.typed"]},
         include_package_data=True,
         data_files=data_files,

symfem/basis_functions.py

@@ -353,8 +353,6 @@ def with_floats(self) -> AnyFunction:
     def __iter__(self) -> typing.Iterator[AnyFunction]:
         """Iterate through components of vector function."""
         f = self.get_function()
-        if not hasattr(f, "__iter__"):
-            raise TypeError(f"'{self.__class__.__name__}' object is not iterable")
         return f.__iter__()
 
 

symfem/elements/dpc.py

@@ -6,7 +6,7 @@
 import sympy
 
 from ..finite_element import CiarletElement
-from ..functionals import DotPointEvaluation, ListOfFunctionals, PointEvaluation
+from ..functionals import DotPointEvaluation, IntegralAgainst, ListOfFunctionals, PointEvaluation
 from ..functions import FunctionInput
 from ..polynomials import polynomial_set_1d, polynomial_set_vector
 from ..references import NonDefaultReferenceError, Reference
@@ -24,20 +24,28 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
             order: The polynomial order
             variant: The variant of the element
         """
+        dofs: ListOfFunctionals = []
         if reference.name == "interval":
-            points = [d.dof_point() for d in Lagrange(reference, order, variant).dofs]
-        elif order == 0:
-            points = [reference.get_point(tuple(
-                sympy.Rational(1, 2) for _ in range(reference.tdim)))]
+            for d in Lagrange(reference, order, variant).dofs:
+                if isinstance(d, PointEvaluation):
+                    dofs.append(PointEvaluation(reference, d.point, entity=(reference.tdim, 0)))
+                elif isinstance(d, IntegralAgainst):
+                    dofs.append(IntegralAgainst(
+                        reference, d.f * reference.jacobian(), entity=(reference.tdim, 0)))
         else:
-            points = [
-                reference.get_point(tuple(sympy.Rational(j, order) for j in i[::-1]))
-                for i in product(range(order + 1), repeat=reference.tdim)
-                if sum(i) <= order
-            ]
+            if order == 0:
+                points = [reference.get_point(tuple(
+                    sympy.Rational(1, 2) for _ in range(reference.tdim)))]
+            else:
+                points = [
+                    reference.get_point(tuple(sympy.Rational(j, order) for j in i[::-1]))
+                    for i in product(range(order + 1), repeat=reference.tdim)
+                    if sum(i) <= order
+                ]
+
+            dofs = [
+                PointEvaluation(reference, d, entity=(reference.tdim, 0)) for d in points]
 
-        dofs: ListOfFunctionals = [
-            PointEvaluation(reference, d, entity=(reference.tdim, 0)) for d in points]
         poly: typing.List[FunctionInput] = []
         poly += polynomial_set_1d(reference.tdim, order)
         poly = reference.map_polyset_from_default(poly)
@@ -59,7 +67,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["interval", "quadrilateral", "hexahedron"]
     min_order = 0
     continuity = "L2"
-    last_updated = "2023.06"
+    last_updated = "2023.07.1"
 
 
 class VectorDPC(CiarletElement):
@@ -107,4 +115,4 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["quadrilateral", "hexahedron"]
     min_order = 0
     continuity = "L2"
-    last_updated = "2023.05"
+    last_updated = "2023.07"

symfem/elements/lagrange.py

@@ -6,9 +6,10 @@
 import sympy
 
 from ..finite_element import CiarletElement
-from ..functionals import DotPointEvaluation, ListOfFunctionals, PointEvaluation
+from ..functionals import DotPointEvaluation, IntegralAgainst, ListOfFunctionals, PointEvaluation
 from ..functions import FunctionInput
-from ..polynomials import polynomial_set_1d, polynomial_set_vector
+from ..polynomials import (lobatto_dual_basis, orthonormal_basis, polynomial_set_1d,
+                           polynomial_set_vector)
 from ..quadrature import get_quadrature
 from ..references import Reference
 
@@ -25,21 +26,34 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
             variant: The variant of the element
         """
         dofs: ListOfFunctionals = []
-        if order == 0:
-            dofs = [
+        if variant == "legendre":
+            basis = orthonormal_basis(reference.name, order, 0)[0]
+            for f in basis:
+                dofs.append(IntegralAgainst(reference, f, (reference.tdim, 0)))
+        elif order == 0:
+            dofs.append(
                 PointEvaluation(
                     reference, reference.get_point(tuple(
                         sympy.Rational(1, reference.tdim + 1)
                         for i in range(reference.tdim)
                     )),
                     entity=(reference.tdim, 0)
                 )
-            ]
+            )
+        elif variant == "lobatto":
+            for v_n, v in enumerate(reference.vertices):
+                dofs.append(PointEvaluation(reference, v, entity=(0, v_n)))
+            for edim in range(1, 4):
+                for e_n in range(reference.sub_entity_count(edim)):
+                    entity = reference.sub_entity(edim, e_n)
+                    basis = lobatto_dual_basis(entity.name, order, False)
+                    for f in basis:
+                        dofs.append(IntegralAgainst(reference, f, (edim, e_n)))
         else:
             points, _ = get_quadrature(variant, order + 1)
 
-            for v_n, v in enumerate(reference.reference_vertices):
-                dofs.append(PointEvaluation(reference, reference.get_point(v), entity=(0, v_n)))
+            for v_n, v in enumerate(reference.vertices):
+                dofs.append(PointEvaluation(reference, v, entity=(0, v_n)))
             for edim in range(1, 4):
                 for e_n in range(reference.sub_entity_count(edim)):
                     entity = reference.sub_entity(edim, e_n)
@@ -66,7 +80,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["interval", "triangle", "tetrahedron"]
     min_order = 0
     continuity = "C0"
-    last_updated = "2023.06"
+    last_updated = "2023.07"
 
 
 class VectorLagrange(CiarletElement):
@@ -85,7 +99,11 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         poly: typing.List[FunctionInput] = []
         if reference.tdim == 1:
             for p in scalar_space.dofs:
-                dofs.append(PointEvaluation(reference, p.dof_point(), entity=p.entity))
+                if isinstance(p, PointEvaluation):
+                    dofs.append(PointEvaluation(reference, p.dof_point(), entity=p.entity))
+                elif isinstance(p, IntegralAgainst):
+                    dofs.append(IntegralAgainst(
+                        reference, p.f * reference.jacobian(), entity=p.entity))
 
             poly += polynomial_set_1d(reference.tdim, order)
         else:
@@ -95,9 +113,15 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
             ]
             for p in scalar_space.dofs:
                 for d in directions:
-                    dofs.append(DotPointEvaluation(reference, p.dof_point(), d, entity=p.entity))
+                    if isinstance(p, PointEvaluation):
+                        dofs.append(DotPointEvaluation(
+                            reference, p.dof_point(), d, entity=p.entity))
+                    elif isinstance(p, IntegralAgainst):
+                        dofs.append(IntegralAgainst(
+                            reference, tuple(p.f * i for i in d), entity=p.entity))
 
             poly += polynomial_set_vector(reference.tdim, reference.tdim, order)
+
         super().__init__(reference, order, poly, dofs, reference.tdim, reference.tdim)
         self.variant = variant
 
@@ -113,7 +137,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["interval", "triangle", "tetrahedron"]
     min_order = 0
     continuity = "C0"
-    last_updated = "2023.05"
+    last_updated = "2023.07.1"
 
 
 class MatrixLagrange(CiarletElement):
@@ -159,7 +183,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["triangle", "tetrahedron"]
     min_order = 0
     continuity = "L2"
-    last_updated = "2023.05"
+    last_updated = "2023.07"
 
 
 class SymmetricMatrixLagrange(CiarletElement):

symfem/elements/lagrange_prism.py

@@ -6,9 +6,9 @@
 import sympy
 
 from ..finite_element import CiarletElement
-from ..functionals import DotPointEvaluation, ListOfFunctionals, PointEvaluation
+from ..functionals import DotPointEvaluation, IntegralAgainst, ListOfFunctionals, PointEvaluation
 from ..functions import FunctionInput
-from ..polynomials import prism_polynomial_set_1d, prism_polynomial_set_vector
+from ..polynomials import orthonormal_basis, prism_polynomial_set_1d, prism_polynomial_set_vector
 from ..quadrature import get_quadrature
 from ..references import NonDefaultReferenceError, Reference
 
@@ -28,7 +28,11 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
             raise NonDefaultReferenceError()
 
         dofs: ListOfFunctionals = []
-        if order == 0:
+        if variant == "legendre":
+            basis = orthonormal_basis(reference.name, order, 0)[0]
+            for f in basis:
+                dofs.append(IntegralAgainst(reference, f, (reference.tdim, 0)))
+        elif order == 0:
             dofs = [
                 PointEvaluation(
                     reference, tuple(
@@ -38,6 +42,8 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
                     entity=(reference.tdim, 0)
                 )
             ]
+        elif variant == "lobatto":
+            raise NotImplementedError()
         else:
             points, _ = get_quadrature(variant, order + 1)
 
@@ -91,7 +97,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["prism"]
     min_order = 0
     continuity = "C0"
-    last_updated = "2023.05"
+    last_updated = "2023.07"
 
 
 class VectorLagrange(CiarletElement):

symfem/elements/lagrange_pyramid.py

@@ -10,9 +10,9 @@
 import sympy
 
 from ..finite_element import CiarletElement
-from ..functionals import ListOfFunctionals, PointEvaluation
+from ..functionals import IntegralAgainst, ListOfFunctionals, PointEvaluation
 from ..functions import FunctionInput
-from ..polynomials import pyramid_polynomial_set_1d
+from ..polynomials import orthonormal_basis, pyramid_polynomial_set_1d
 from ..quadrature import get_quadrature
 from ..references import NonDefaultReferenceError, Reference
 
@@ -32,7 +32,11 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
             raise NonDefaultReferenceError()
 
         dofs: ListOfFunctionals = []
-        if order == 0:
+        if variant == "legendre":
+            basis = orthonormal_basis(reference.name, order, 0)[0]
+            for f in basis:
+                dofs.append(IntegralAgainst(reference, f, (reference.tdim, 0)))
+        elif order == 0:
             dofs = [
                 PointEvaluation(
                     reference, tuple(
@@ -42,6 +46,8 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
                     entity=(reference.tdim, 0)
                 )
             ]
+        elif variant == "lobatto":
+            raise NotImplementedError()
         else:
             points, _ = get_quadrature(variant, order + 1)
 
@@ -97,4 +103,4 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["pyramid"]
     min_order = 0
     continuity = "C0"
-    last_updated = "2023.05"
+    last_updated = "2023.07"

symfem/elements/q.py

@@ -6,12 +6,12 @@
 import sympy
 
 from ..finite_element import CiarletElement, FiniteElement
-from ..functionals import (DotPointEvaluation, IntegralMoment, ListOfFunctionals,
+from ..functionals import (DotPointEvaluation, IntegralAgainst, IntegralMoment, ListOfFunctionals,
                            NormalIntegralMoment, PointEvaluation, TangentIntegralMoment)
 from ..functions import FunctionInput
 from ..moments import make_integral_moment_dofs
-from ..polynomials import (Hcurl_quolynomials, Hdiv_quolynomials, quolynomial_set_1d,
-                           quolynomial_set_vector)
+from ..polynomials import (Hcurl_quolynomials, Hdiv_quolynomials, lobatto_dual_basis,
+                           orthonormal_basis, quolynomial_set_1d, quolynomial_set_vector)
 from ..quadrature import get_quadrature
 from ..references import NonDefaultReferenceError, Reference
 
@@ -28,11 +28,26 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
             variant: The variant of the element
         """
         dofs: ListOfFunctionals = []
-        if order == 0:
+        if variant == "legendre":
+            basis = orthonormal_basis(reference.name, order, 0)[0]
+            for f in basis:
+                dofs.append(IntegralAgainst(reference, f, (reference.tdim, 0)))
+        elif order == 0:
             dofs = [PointEvaluation(
                 reference,
                 reference.get_point(tuple(sympy.Rational(1, 2) for i in range(reference.tdim))),
                 entity=(reference.tdim, 0))]
+        elif variant == "lobatto":
+            if reference != reference.default_reference():
+                raise NonDefaultReferenceError()
+            for v_n, v in enumerate(reference.vertices):
+                dofs.append(PointEvaluation(reference, v, entity=(0, v_n)))
+            for edim in range(1, 4):
+                for e_n in range(reference.sub_entity_count(edim)):
+                    entity = reference.sub_entity(edim, e_n)
+                    basis = lobatto_dual_basis(entity.name, order, False)
+                    for f in basis:
+                        dofs.append(IntegralAgainst(reference, f, (edim, e_n)))
         else:
             points, _ = get_quadrature(variant, order + 1)
 
@@ -64,6 +79,8 @@ def get_tensor_factorisation(
         Returns:
             The tensor factorisation
         """
+        if self.variant == "lobatto":
+            return super().get_tensor_factorisation()
         from symfem import create_element
         interval_q = create_element("interval", "Lagrange", self.order)
 
@@ -110,7 +127,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["quadrilateral", "hexahedron"]
     min_order = 0
     continuity = "C0"
-    last_updated = "2023.06"
+    last_updated = "2023.07.1"
 
 
 class VectorQ(CiarletElement):

symfem/elements/serendipity.py

@@ -103,7 +103,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["quadrilateral", "hexahedron"]
     min_order = 1
     continuity = "H(curl)"
-    last_updated = "2023.06"
+    last_updated = "2023.07"
 
 
 class SerendipityDiv(CiarletElement):
@@ -146,4 +146,4 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["quadrilateral", "hexahedron"]
     min_order = 1
     continuity = "H(div)"
-    last_updated = "2023.06.1"
+    last_updated = "2023.07"

symfem/functions.py

@@ -340,6 +340,10 @@ def with_floats(self) -> AnyFunction:
         """
         pass
 
+    def __iter__(self) -> typing.Iterator[AnyFunction]:
+        """Iterate through components of vector function."""
+        raise TypeError(f"'{self.__class__.__name__}' object is not iterable")
+
     def integrate(self, *limits: typing.Tuple[
         sympy.core.symbol.Symbol,
         typing.Union[int, sympy.core.expr.Expr],