Update HHJ construction on tetrahedra and fix map (#267)

* Update tet construction and fix map * Fix typo * Fix typo * Fix parenthesis
mscroggs · May 7, 2024 · cb07771 · cb07771
1 parent eb9b1e7
commit cb07771
Showing 1 changed file with 9 additions and 7 deletions.
diff --git a/symfem/elements/hhj.py b/symfem/elements/hhj.py
@@ -5,6 +5,8 @@
 
 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
+or (Pechstein, Schöberl, 2018) for a more recent version for tetrahedra
+https://doi.org/10.1007/s00211-017-0933-3
 """
 
 import typing
@@ -47,10 +49,10 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         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 = [((0, 1, 1), (1, 0, 1), (1, 1, 0)),
+                          ((-6, 1, 1), (1, 0, 1), (1, 1, 0)),
+                          ((0, 1, 1), (1, -6, 1), (1, 1, 0)),
+                          ((0, 1, 1), (1, 0, 1), (1, 1, -6))]
             directions_extra = [((0, 0, -1), (0, 0, 1), (-1, 1, 0)),
                                 ((0, -1, 0), (-1, 0, 1), (0, 1, 0))]
 
@@ -68,15 +70,15 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
                 for d in directions:
                     dofs.append(IntegralMoment(
                         reference, p * MatrixFunction(d), dof, entity=(reference.tdim, 0),
-                        mapping="double_covariant"))
+                        mapping="double_contravariant"))
         # 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"))
+                    mapping="double_contravariant"))
 
         self.variant = variant
 
@@ -95,4 +97,4 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     references = ["triangle", "tetrahedron"]
     min_order = 0
     continuity = "inner H(div)"
-    last_updated = "2023.08"
+    last_updated = "2024.02"