Add Lederer-Schoberl element (#301)

* add lowest degree Lederer-Schoberl element

* ruff

* higher order element on triangles

* Fix element in 3D

* mypy

* update docs

* co_contravariant

* cahce

* changelog

Author: mscroggs

CHANGELOG_SINCE_LAST_VERSION.md

@@ -1,3 +1,4 @@
 - Allow "gl" variant of Lagrange element
 - Implement component-wise integrals of VectorFunction and MatrixFunction
 - Implement grad of VectorFunction
+- Add Lederer-Schoberl element

README.md

@@ -187,6 +187,7 @@ The reference triangle has vertices (0, 0), (1, 0), and (0, 1). Its sub-entities
 - Hsieh-Clough-Tocher (alternative names: Clough-Tocher, HCT, CT)
 - Kong-Mulder-Veldhuizen (alternative names: KMV)
 - Lagrange (alternative names: P)
+- Lederer-Schoberl
 - Mardal-Tai-Winther (alternative names: MTW)
 - matrix Lagrange
 - Morley
@@ -258,6 +259,7 @@ The reference tetrahedron has vertices (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0,
 - Hermite
 - Kong-Mulder-Veldhuizen (alternative names: KMV)
 - Lagrange (alternative names: P)
+- Lederer-Schoberl
 - Mardal-Tai-Winther (alternative names: MTW)
 - matrix Lagrange
 - Morley-Wang-Xu (alternative names: MWX)

symfem/create.py

@@ -178,7 +178,8 @@ def create_element(
                       enriched Galerkin, EG,
                       enriched vector Galerkin, locking-free enriched Galerkin, LFEG,
                       P1 macro,
-                      Alfeld-Sorokina, AS
+                      Alfeld-Sorokina, AS,
+                      Lederer-Schoberl
         order: The order of the element.
         vertices: The vertices of the reference.
     """

symfem/elements/lederer_schoberl.py

@@ -0,0 +1,164 @@
+"""Lederer-Schoberl elements on simplices.
+
+This element's definition appears in https://doi.org/10.34726/hss.2019.62042
+(Lederer, 2019).
+"""
+
+import typing
+
+
+from symfem.elements.lagrange import Lagrange
+from symfem.finite_element import CiarletElement
+from symfem.functionals import (
+    InnerProductIntegralMoment,
+    IntegralMoment,
+    TraceIntegralMoment,
+    ListOfFunctionals,
+)
+from symfem.functions import FunctionInput, MatrixFunction
+from symfem.moments import make_integral_moment_dofs
+from symfem.polynomials import polynomial_set_vector
+from symfem.references import Reference
+from symfem.symbols import t, x
+
+__all__ = ["LedererSchoberl"]
+
+
+class LedererSchoberl(CiarletElement):
+    """A Lederer-Schoberl element on a simplex."""
+
+    def __init__(self, reference: Reference, order: int):
+        """Create the element.
+
+        Args:
+            reference: The reference element
+            order: The polynomial order
+            variant: The variant of the element
+        """
+        from symfem import create_reference
+
+        tdim = reference.tdim
+        poly: typing.List[FunctionInput] = [
+            tuple(tuple(p[i * tdim : (i + 1) * tdim]) for i in range(tdim))
+            for p in polynomial_set_vector(tdim, tdim**2, order)
+        ]
+
+        dofs: ListOfFunctionals = []
+        facet_dim = reference.tdim - 1
+        space = Lagrange(
+            create_reference(["point", "interval", "triangle"][facet_dim]), order, "equispaced"
+        )
+        basis = [f.subs(x, t) for f in space.get_basis_functions()]
+        for facet_n in range(reference.sub_entity_count(facet_dim)):
+            facet = reference.sub_entity(facet_dim, facet_n)
+            for tangent in facet.axes:
+                for f, dof in zip(basis, space.dofs):
+                    dofs.append(
+                        InnerProductIntegralMoment(
+                            reference,
+                            f,
+                            tuple(i / facet.volume() for i in tangent),
+                            facet.normal(),
+                            dof,
+                            entity=(facet_dim, facet_n),
+                            mapping="co_contravariant",
+                        )
+                    )
+
+        dofs += make_integral_moment_dofs(
+            reference,
+            cells=(
+                TraceIntegralMoment,
+                Lagrange,
+                order,
+                "co_contravariant",
+            ),
+        )
+
+        if order > 0:
+            if reference.tdim == 2:
+                functions = [
+                    MatrixFunction(i)
+                    for i in [
+                        (((x[0] + x[1] - 1) / 2, 0), (0, (1 - x[0] - x[1]) / 2)),
+                        ((x[0] / 2, 0), (x[0], -x[0] / 2)),
+                        ((x[1] / 2, -x[1]), (0, -x[1] / 2)),
+                    ]
+                ]
+            else:
+                assert reference.tdim == 3
+                functions = [
+                    MatrixFunction(i)
+                    for i in [
+                        (
+                            (2 * (1 - x[0] - x[1] - x[2]) / 3, 0, 0),
+                            (0, (x[0] + x[1] + x[2] - 1) / 3, 0),
+                            (0, 0, (x[0] + x[1] + x[2] - 1) / 3),
+                        ),
+                        (
+                            ((x[0] + x[1] + x[2] - 1) / 3, 0, 0),
+                            (0, 2 * (1 - x[0] - x[1] - x[2]) / 3, 0),
+                            (0, 0, (x[0] + x[1] + x[2] - 1) / 3),
+                        ),
+                        ((x[0] / 3, 0, 0), (x[0], -2 * x[0] / 3, 0), (0, 0, x[0] / 3)),
+                        ((x[0] / 3, 0, 0), (0, x[0] / 3, 0), (x[0], 0, -2 * x[0] / 3)),
+                        ((-x[1] / 3, 0, 0), (0, -x[1] / 3, 0), (0, -x[1], 2 * x[1] / 3)),
+                        ((-x[1] / 3, x[1], 0), (0, 2 * x[1] / 3, 0), (0, x[1], -x[1] / 3)),
+                        ((x[2] / 3, 0, -x[2]), (0, x[2] / 3, -x[2]), (0, 0, -2 * x[2] / 3)),
+                        ((-2 * x[2] / 3, 0, x[2]), (0, x[2] / 3, 0), (0, 0, x[2] / 3)),
+                    ]
+                ]
+
+            lagrange = Lagrange(reference, order - 1)
+            for dof, f in zip(lagrange.dofs, lagrange.get_basis_functions()):
+                for g in functions:
+                    dofs.append(
+                        IntegralMoment(
+                            reference,
+                            f * g,
+                            dof,
+                            (reference.tdim, 0),
+                            "co_contravariant",
+                        )
+                    )
+
+        super().__init__(
+            reference,
+            order,
+            poly,
+            dofs,
+            reference.tdim,
+            reference.tdim**2,
+            (reference.tdim, reference.tdim),
+        )
+
+    def init_kwargs(self) -> typing.Dict[str, typing.Any]:
+        """Return the kwargs used to create this element.
+
+        Returns:
+            Keyword argument dictionary
+        """
+        return {}
+
+    @property
+    def lagrange_subdegree(self) -> int:
+        return self.order
+
+    @property
+    def lagrange_superdegree(self) -> typing.Optional[int]:
+        return self.order
+
+    @property
+    def polynomial_subdegree(self) -> int:
+        return self.order
+
+    @property
+    def polynomial_superdegree(self) -> typing.Optional[int]:
+        return self.order
+
+    names = ["Lederer-Schoberl"]
+    references = ["triangle", "tetrahedron"]
+    min_order = 0
+    continuity = "inner H(curl div)"
+    value_type = "matrix"
+    last_updated = "2024.03"

symfem/finite_element.py

@@ -521,6 +521,13 @@ def get_piece(f, point):
                             assert dim == 2
                             f = [f[1, 1], f[2, 2]]
                             g = [g[1, 1], g[2, 2]]
+                elif continuity == "inner H(curl div)":
+                    if len(vertices[0]) == 2:
+                        f = f[1, 0]
+                        g = g[1, 0]
+                    if len(vertices[0]) == 3:
+                        f = [f[1, 0], f[2, 0]]
+                        g = [g[1, 0], g[2, 0]]
                 elif continuity == "integral inner H(div)":
                     f = f[0, 0].integrate((x[1], 0, 1))
                     g = g[0, 0].integrate((x[1], 0, 1))

symfem/functionals.py

@@ -39,6 +39,7 @@
     "DivergenceIntegralMoment",
     "TangentIntegralMoment",
     "NormalIntegralMoment",
+    "TraceIntegralMoment",
     "NormalDerivativeIntegralMoment",
     "InnerProductIntegralMoment",
     "NormalInnerProductIntegralMoment",
@@ -1638,6 +1639,68 @@ def get_tex(self) -> typing.Tuple[str, typing.List[str]]:
     name = "Normal integral moment"
 
 
+class TraceIntegralMoment(IntegralMoment):
+    """An integral moment of the trace of a matrix."""
+
+    def __init__(
+        self,
+        reference: Reference,
+        f_in: FunctionInput,
+        dof: BaseFunctional,
+        entity: typing.Tuple[int, int],
+        mapping: typing.Union[str, None] = "contravariant",
+    ):
+        """Create the functional.
+
+        Args:
+            reference: The reference cell
+            f_in: The function to multiply with
+            dof: The DOF in a moment space that the function is associated with
+            entity: The entity the functional is associated with
+            mapping: The type of mappping from the reference cell to a physical cell
+        """
+        scalar_f = parse_function_input(f_in).subs(x, t)
+        assert scalar_f.is_scalar
+        super().__init__(
+            reference,
+            scalar_f,
+            dof,
+            entity=entity,
+            mapping=mapping,
+            map_function=False,
+        )
+
+    def dot(self, function: Function) -> ScalarFunction:
+        """Take the product of a function with the trace of the matrix.
+
+        Args:
+            function: The function
+
+        Returns:
+            The inner product of the function and the matrix's trace
+        """
+        assert function.is_matrix and function.shape[0] == function.shape[1]
+        trace = sum(function[i, i] for i in range(function.shape[0]))
+        return trace * self.f * self.integral_domain.jacobian()
+
+    def get_tex(self) -> typing.Tuple[str, typing.List[str]]:
+        """Get a representation of the functional as TeX, and list of terms involved.
+
+        Returns:
+            Representation of the functional as TeX, and list of terms involved
+        """
+        entity = self.entity_tex()
+        entity_def = self.entity_definition()
+        desc = "\\boldsymbol{V}\\mapsto"
+        desc += f"\\displaystyle\\int_{{{entity}}}"
+        desc += "\\operatorname{tr}(\\boldsymbol{V})"
+        if self.f != 1:
+            desc += "(" + _to_tex(self.f, True) + ")"
+        return desc, [entity_def]
+
+    name = "Trace integral moment"
+
+
 class NormalDerivativeIntegralMoment(DerivativeIntegralMoment):
     """An integral moment in the normal direction."""