Add Sinwel construction for HHJ

symfem/elements/hhj.py

@@ -2,17 +2,20 @@
 
 This element's definition appears in https://arxiv.org/abs/1909.09687
 (Arnold, Walker, 2020)
+
+For an alternative construction see (Sinwel, 2009) and sections 4.4.2.2 and 4.4.3.2
+https://numa.jku.at/media/filer_public/b7/42/b74263c9-f723-4076-b1b2-c2726126bf32/phd-sinwel.pdf
 """
 
 import typing
 
 from ..finite_element import CiarletElement
 from ..functionals import IntegralMoment, ListOfFunctionals, NormalInnerProductIntegralMoment
-from ..functions import FunctionInput
+from ..functions import FunctionInput, MatrixFunction
 from ..moments import make_integral_moment_dofs
 from ..polynomials import polynomial_set_vector
 from ..references import NonDefaultReferenceError, Reference
-from .lagrange import Lagrange, SymmetricMatrixLagrange
+from .lagrange import Lagrange
 
 
 class HellanHerrmannJohnson(CiarletElement):
@@ -28,20 +31,53 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         """
         if reference.vertices != reference.reference_vertices:
             raise NonDefaultReferenceError()
-        assert reference.name == "triangle"
+        assert reference.name == "triangle" or reference.name == "tetrahedron"
 
         poly: typing.List[FunctionInput] = []
-        poly += [((p[0], p[1]), (p[1], p[2]))
-                 for p in polynomial_set_vector(reference.tdim, 3, order)]
+        directions: typing.List[typing.Tuple[typing.Tuple[int, ...], ...]] = []
+        directions_extra: typing.List[typing.Tuple[typing.Tuple[int, ...], ...]] = []
+
+        if reference.tdim == 2:
+            poly = [((p[0], p[1]), (p[1], p[2]))
+                    for p in polynomial_set_vector(reference.tdim, 3, order)]
+            directions = [((0, 1), (1, 0)),
+                          ((-2, 1), (1, 0)),
+                          ((0, -1), (-1, 2))]
+            directions_extra = []
+        if reference.tdim == 3:
+            poly = [((p[0], p[1], p[2]), (p[1], p[3], p[4]), (p[2], p[4], p[5]))
+                    for p in polynomial_set_vector(reference.tdim, 6, order)]
+            directions = [((-2, 1, 0), (1, 0, 0), (0, 0, 0)),
+                          ((0, 1, -1), (1, -2, 1), (-1, 1, 0)),
+                          ((0, 0, 0), (0, 0, -1), (0, -1, 2)),
+                          ((0, 0, 1), (0, 0, 0), (1, 0, 0))]
+            directions_extra = [((0, 0, -1), (0, 0, 1), (-1, 1, 0)),
+                                ((0, -1, 0), (-1, 0, 1), (0, 1, 0))]
 
         dofs: ListOfFunctionals = make_integral_moment_dofs(
             reference,
             facets=(NormalInnerProductIntegralMoment, Lagrange, order,
                     {"variant": variant}),
-            cells=(IntegralMoment, SymmetricMatrixLagrange, order - 1,
-                   {"variant": variant}),
         )
 
+        # cell functions
+        if order > 0:
+            space = Lagrange(reference, order - 1, variant)
+            basis = space.get_basis_functions()
+            for p, dof in zip(basis, space.dofs):
+                for d in directions:
+                    dofs.append(IntegralMoment(
+                        reference, p * MatrixFunction(d), dof, entity=(reference.tdim, 0),
+                        mapping="double_covariant"))
+        # cell functions extra
+        space_extra = Lagrange(reference, order, variant)
+        basis_extra = space_extra.get_basis_functions()
+        for p, dof in zip(basis_extra, space_extra.dofs):
+            for d in directions_extra:
+                dofs.append(IntegralMoment(
+                    reference, p * MatrixFunction(d), dof, entity=(reference.tdim, 0),
+                    mapping="double_covariant"))
+
         self.variant = variant
 
         super().__init__(reference, order, poly, dofs, reference.tdim, reference.tdim ** 2,
@@ -56,7 +92,7 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
         return {"variant": self.variant}
 
     names = ["Hellan-Herrmann-Johnson", "HHJ"]
-    references = ["triangle"]
+    references = ["triangle", "tetrahedron"]
     min_order = 0
     continuity = "inner H(div)"
-    last_updated = "2023.06"
+    last_updated = "2023.08"

test/test_hellan_herrmann_johnson.py

@@ -0,0 +1,41 @@
+"""Test Hellan-Herrmann-Johnson element."""
+
+import pytest
+
+from symfem import create_element
+from symfem.symbols import x
+from symfem.utils import allequal
+from symfem.functions import VectorFunction
+
+
+@pytest.mark.parametrize('reference', ['triangle', 'tetrahedron'])
+@pytest.mark.parametrize('order', [0, 1, 2])
+def test_create(reference, order):
+
+    element = create_element(reference, "HHJ", order)
+
+    # Get the basis functions associated with the interior
+    basis = element.get_basis_functions()
+    functions = [basis[d] for d in element.entity_dofs(element.reference.tdim, 0)]
+
+    if reference == "triangle":
+        # Check that these tensor functions have zero normal normal component on each edge
+        for f in functions:
+            M, n = f.subs(x[0], 1 - x[1]), VectorFunction((1, 1))
+            assert allequal((M @ n).dot(n), 0)
+            M, n = f.subs(x[0], 0), VectorFunction((-1, 0))
+            assert allequal((M @ n).dot(n), 0)
+            M, n = f.subs(x[1], 0), VectorFunction((0, -1))
+            assert allequal((M @ n).dot(n), 0)
+
+    if reference == "tetrahedron":
+        # Check that these tensor functions have zero normal normal component on each face
+        for f in functions:
+            M, n = f.subs(x[0], 1 - x[1] - x[2]), VectorFunction((1, 1, 1))
+            assert allequal((M @ n).dot(n), 0)
+            M, n = f.subs(x[0], 0), VectorFunction((-1, 0, 0))
+            assert allequal((M @ n).dot(n), 0)
+            M, n = f.subs(x[1], 0), VectorFunction((0, -1, 0))
+            assert allequal((M @ n).dot(n), 0)
+            M, n = f.subs(x[2], 0), VectorFunction((0, 0, -1))
+            assert allequal((M @ n).dot(n), 0)