Merge branch 'main' into mscroggs/bernardi-raugel

CHANGELOG_SINCE_LAST_VERSION.md

@@ -1,2 +1,4 @@
 - Corrected nonconforming Arnold-Winther polyset
 - Correct Alfeld-Sorokina range dim
+- Correct Guzman-Neilan second kind element
+- Add Guzman-Neilan first kind element

README.md

@@ -180,7 +180,8 @@ The reference triangle has vertices (0, 0), (1, 0), and (0, 1). Its sub-entities
 - enriched Galerkin (alternative names: EG)
 - enriched vector Galerkin (alternative names: locking-free enriched Galerkin, LFEG)
 - Fortin-Soulie (alternative names: FS)
-- Guzman-Neilan
+- Guzman-Neilan first kind (alternative names: Guzman-Neilan)
+- Guzman-Neilan second kind
 - Hellan-Herrmann-Johnson (alternative names: HHJ)
 - Hermite
 - Hsieh-Clough-Tocher (alternative names: Clough-Tocher, HCT, CT)
@@ -251,7 +252,8 @@ The reference tetrahedron has vertices (0, 0, 0), (1, 0, 0), (0, 1, 0), and (0,
 - Crouzeix-Raviart (alternative names: CR, Crouzeix-Falk, CF)
 - enriched Galerkin (alternative names: EG)
 - enriched vector Galerkin (alternative names: locking-free enriched Galerkin, LFEG)
-- Guzman-Neilan
+- Guzman-Neilan first kind (alternative names: Guzman-Neilan)
+- Guzman-Neilan second kind
 - Hellan-Herrmann-Johnson (alternative names: HHJ)
 - Hermite
 - Kong-Mulder-Veldhuizen (alternative names: KMV)

symfem/create.py

@@ -161,7 +161,8 @@ def create_element(
                       Bernardi-Raugel,
                       Wu-Xu,
                       transition,
-                      Guzman-Neilan,
+                      Guzman-Neilan first kind, Guzman-Neilan,
+                      Guzman-Neilan second kind,
                       nonconforming Arnold-Winther, nonconforming AW,
                       TScurl, trimmed serendipity Hcurl,
                       TSdiv, trimmed serendipity Hdiv,

symfem/elements/guzman_neilan.py

@@ -18,11 +18,157 @@
 from symfem.elements.bernardi_raugel import BernardiRaugel
 from symfem.elements.lagrange import Lagrange, VectorLagrange
 
-__all__ = ["GuzmanNeilan", "make_piecewise_lagrange"]
+__all__ = ["GuzmanNeilanFirstKind", "GuzmanNeilanSecondKind", "make_piecewise_lagrange"]
 
 
-class GuzmanNeilan(CiarletElement):
-    """Guzman-Neilan Hdiv finite element."""
+class GuzmanNeilanFirstKind(CiarletElement):
+    """Guzman-Neilan finite element."""
+
+    def __init__(self, reference: Reference, order: int):
+        """Create the element.
+
+        Args:
+            reference: The reference element
+            order: The polynomial order
+        """
+        if reference.vertices != reference.reference_vertices:
+            raise NonDefaultReferenceError()
+
+        if reference.name == "triangle":
+            poly = self._make_polyset_triangle(reference, order)
+        else:
+            poly = self._make_polyset_tetrahedron(reference, order)
+
+        br = BernardiRaugel(reference, order)
+        if order == 1:
+            dofs = br.dofs
+        else:
+            assert order == 2 and reference.name == "tetrahedron"
+            dofs = br.dofs[:-3]
+
+        super().__init__(reference, order, poly, dofs, reference.tdim, reference.tdim)
+
+    def _make_polyset_triangle(
+        self, reference: Reference, order: int
+    ) -> typing.List[FunctionInput]:
+        """Make the polyset for a triangle.
+
+        Args:
+            reference: The reference cell
+            order: The polynomial order
+
+        Returns:
+            The polynomial set
+        """
+        assert order == 1
+
+        from symfem.elements._guzman_neilan_triangle import coeffs
+
+        mid = reference.midpoint()
+
+        sub_tris: typing.List[SetOfPoints] = [
+            (reference.vertices[0], reference.vertices[1], mid),
+            (reference.vertices[0], reference.vertices[2], mid),
+            (reference.vertices[1], reference.vertices[2], mid),
+        ]
+
+        lagrange = VectorLagrange(reference, order)
+        basis: typing.List[FunctionInput] = [
+            PiecewiseFunction({i: p for i in sub_tris}, 2) for p in lagrange.get_polynomial_basis()
+        ]
+
+        sub_basis = make_piecewise_lagrange(sub_tris, "triangle", reference.tdim, True)
+        fs = BernardiRaugel(reference, 1).get_basis_functions()[-3:]
+        for c, f in zip(coeffs, fs):
+            assert isinstance(f, VectorFunction)
+            pieces: typing.Dict[SetOfPointsInput, FunctionInput] = {}
+            for tri in sub_tris:
+                function: typing.List[sympy.core.expr.Expr] = []
+                for j in range(reference.tdim):
+                    component = f[j]
+                    for k, b in zip(c, sub_basis):
+                        component -= k * b.pieces[tri][j]
+                    c_sym = component.as_sympy()
+                    assert isinstance(c_sym, sympy.core.expr.Expr)
+                    function.append(c_sym)
+                pieces[tri] = VectorFunction(tuple(function))
+            basis.append(PiecewiseFunction(pieces, 2))
+        return basis
+
+    def _make_polyset_tetrahedron(
+        self, reference: Reference, order: int
+    ) -> typing.List[FunctionInput]:
+        """Make the polyset for a tetrahedron.
+
+        Args:
+            reference: The reference cell
+            order: The polynomial order
+
+        Returns:
+            The polynomial set
+        """
+        assert order in [1, 2]
+        from symfem.elements._guzman_neilan_tetrahedron import coeffs
+
+        mid = reference.midpoint()
+
+        sub_tets: typing.List[SetOfPoints] = [
+            (reference.vertices[0], reference.vertices[1], reference.vertices[2], mid),
+            (reference.vertices[0], reference.vertices[1], reference.vertices[3], mid),
+            (reference.vertices[0], reference.vertices[2], reference.vertices[3], mid),
+            (reference.vertices[1], reference.vertices[2], reference.vertices[3], mid),
+        ]
+
+        lagrange = VectorLagrange(reference, order)
+        basis: typing.List[FunctionInput] = [
+            PiecewiseFunction({i: p for i in sub_tets}, 3) for p in lagrange.get_polynomial_basis()
+        ]
+
+        sub_basis = make_piecewise_lagrange(sub_tets, "tetrahedron", reference.tdim, True)
+        fs = BernardiRaugel(reference, 1).get_basis_functions()[-4:]
+        for c, f in zip(coeffs, fs):
+            assert isinstance(f, VectorFunction)
+            pieces: typing.Dict[SetOfPointsInput, FunctionInput] = {}
+            for tet in sub_tets:
+                function: typing.List[sympy.core.expr.Expr] = []
+                for j in range(reference.tdim):
+                    component = f[j]
+                    for k, b in zip(c, sub_basis):
+                        component -= k * b.pieces[tet][j]
+                    c_sym = component.as_sympy()
+                    assert isinstance(c_sym, sympy.core.expr.Expr)
+                    function.append(c_sym)
+                pieces[tet] = VectorFunction(tuple(function))
+            basis.append(PiecewiseFunction(pieces, 3))
+        return basis
+
+    @property
+    def lagrange_subdegree(self) -> int:
+        raise NotImplementedError()
+
+    @property
+    def lagrange_superdegree(self) -> typing.Optional[int]:
+        raise NotImplementedError()
+
+    @property
+    def polynomial_subdegree(self) -> int:
+        raise NotImplementedError()
+
+    @property
+    def polynomial_superdegree(self) -> typing.Optional[int]:
+        raise NotImplementedError()
+
+    names = ["Guzman-Neilan first kind", "Guzman-Neilan"]
+    references = ["triangle", "tetrahedron"]
+    min_order = 1
+    max_order = {"triangle": 1, "tetrahedron": 2}
+    continuity = "L2"
+    value_type = "vector macro"
+    last_updated = "2024.10"
+
+
+class GuzmanNeilanSecondKind(CiarletElement):
+    """Guzman-Neilan second kind finite element."""
 
     def __init__(self, reference: Reference, order: int):
         """Create the element.
@@ -41,16 +187,21 @@ def __init__(self, reference: Reference, order: int):
 
         dofs: ListOfFunctionals = []
 
-        for n in range(reference.sub_entity_count(reference.tdim - 1)):
-            facet = reference.sub_entity(reference.tdim - 1, n)
-            for v in facet.vertices:
+        tdim = reference.tdim
+
+        # Evaluation at vertices
+        for n in range(reference.sub_entity_count(0)):
+            vertex = reference.sub_entity(0, n)
+            v = vertex.vertices[0]
+            for i in range(tdim):
+                direction = tuple(1 if i == j else 0 for j in range(tdim))
                 dofs.append(
                     DotPointEvaluation(
                         reference,
                         v,
-                        tuple(i * facet.jacobian() for i in facet.normal()),
-                        entity=(reference.tdim - 1, n),
-                        mapping="contravariant",
+                        direction,
+                        entity=(0, n),
+                        mapping="identity",
                     )
                 )
 
@@ -60,28 +211,15 @@ def __init__(self, reference: Reference, order: int):
             # Midpoints of edges
             for n in range(reference.sub_entity_count(1)):
                 edge = reference.sub_entity(1, n)
-                dofs.append(
-                    DotPointEvaluation(
-                        reference,
-                        edge.midpoint(),
-                        tuple(i * edge.jacobian() for i in edge.tangent()),
-                        entity=(1, n),
-                        mapping="contravariant",
-                    )
-                )
-
-            # Midpoints of edges of faces
-            for n in range(reference.sub_entity_count(2)):
-                face = reference.sub_entity(2, n)
-                for m in range(3):
-                    edge = face.sub_entity(1, m)
+                for i in range(tdim):
+                    direction = tuple(1 if i == j else 0 for j in range(tdim))
                     dofs.append(
                         DotPointEvaluation(
                             reference,
                             edge.midpoint(),
-                            tuple(i * face.jacobian() for i in face.normal()),
-                            entity=(2, n),
-                            mapping="contravariant",
+                            direction,
+                            entity=(1, n),
+                            mapping="identity",
                         )
                     )
 
@@ -90,7 +228,7 @@ def __init__(self, reference: Reference, order: int):
                 p = tuple((i + sympy.Rational(1, 4)) / 2 for i in v)
                 for d in [(1, 0, 0), (0, 1, 0), (0, 0, 1)]:
                     dofs.append(
-                        DotPointEvaluation(reference, p, d, entity=(3, 0), mapping="contravariant")
+                        DotPointEvaluation(reference, p, d, entity=(3, 0), mapping="identity")
                     )
 
         dofs += make_integral_moment_dofs(
@@ -99,15 +237,13 @@ def __init__(self, reference: Reference, order: int):
         )
 
         mid = reference.midpoint()
-        for i in range(reference.tdim):
-            direction = tuple(1 if i == j else 0 for j in range(reference.tdim))
+        for i in range(tdim):
+            direction = tuple(1 if i == j else 0 for j in range(tdim))
             dofs.append(
-                DotPointEvaluation(
-                    reference, mid, direction, entity=(reference.tdim, 0), mapping="contravariant"
-                )
+                DotPointEvaluation(reference, mid, direction, entity=(tdim, 0), mapping="identity")
             )
 
-        super().__init__(reference, order, poly, dofs, reference.tdim, reference.tdim)
+        super().__init__(reference, order, poly, dofs, tdim, tdim)
 
     def _make_polyset_triangle(
         self, reference: Reference, order: int
@@ -215,13 +351,13 @@ def polynomial_subdegree(self) -> int:
     def polynomial_superdegree(self) -> typing.Optional[int]:
         raise NotImplementedError()
 
-    names = ["Guzman-Neilan"]
+    names = ["Guzman-Neilan second kind"]
     references = ["triangle", "tetrahedron"]
     min_order = 1
     max_order = {"triangle": 1, "tetrahedron": 2}
-    continuity = "H(div)"
+    continuity = "L2"
     value_type = "vector macro"
-    last_updated = "2023.06"
+    last_updated = "2024.10"
 
 
 def make_piecewise_lagrange(

test/test_guzman_neilan.py

@@ -12,7 +12,7 @@
 
 @pytest.mark.parametrize("order", [1])
 def test_guzman_neilan_triangle(order):
-    e = symfem.create_element("triangle", "Guzman-Neilan", order)
+    e = symfem.create_element("triangle", "Guzman-Neilan second kind", order)
 
     for p in e._basis[-3:]:
         for piece in p.pieces.values():
@@ -21,7 +21,7 @@ def test_guzman_neilan_triangle(order):
 
 @pytest.mark.parametrize("order", [1, 2])
 def test_guzman_neilan_tetrahedron(order):
-    e = symfem.create_element("tetrahedron", "Guzman-Neilan", order)
+    e = symfem.create_element("tetrahedron", "Guzman-Neilan second kind", order)
 
     mid = tuple(sympy.Rational(sum(i), len(i)) for i in zip(*e.reference.vertices))
     for p in e._basis[-4:]:
@@ -34,7 +34,7 @@ def test_guzman_neilan_tetrahedron(order):
 @pytest.mark.parametrize("order", [1])
 def test_basis_continuity_triangle(order):
     N = 5
-    e = symfem.create_element("triangle", "Guzman-Neilan", order)
+    e = symfem.create_element("triangle", "Guzman-Neilan second kind", order)
     third = sympy.Rational(1, 3)
     one = sympy.Integer(1)
     for pt in [(0, 0), (1, 0), (0, 1), (third, third)]:
@@ -60,7 +60,7 @@ def test_basis_continuity_triangle(order):
 @pytest.mark.parametrize("order", [1, 2])
 def test_basis_continuity_tetrahedron(order):
     N = 5
-    e = symfem.create_element("tetrahedron", "Guzman-Neilan", order)
+    e = symfem.create_element("tetrahedron", "Guzman-Neilan second kind", order)
     quarter = sympy.Rational(1, 4)
     one = sympy.Integer(1)
     for pt in [(0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1), (quarter, quarter, quarter)]:

test/test_plotting.py

@@ -204,7 +204,7 @@ def test_function_plots_piecewise_scalar(reference):
     ],
 )
 def test_function_plots_piecewise_vector(reference):
-    e = symfem.create_element(reference, "Guzman-Neilan", 1)
+    e = symfem.create_element(reference, "Guzman-Neilan second kind", 1)
     for ext in ["svg", "png", "tex"]:
         e.plot_basis_function(
             0, os.path.join(folder, f"test_function_plots_piecewise_vector-{reference}.{ext}")

test/utils.py

@@ -63,7 +63,8 @@
         "FS": [({}, [2])],
         "Taylor": [({}, range(0, 5))],
         "Bernardi-Raugel": [({}, [1])],
-        "Guzman-Neilan": [({}, [1])],
+        "Guzman-Neilan first kind": [({}, [1])],
+        "Guzman-Neilan second kind": [({}, [1])],
         "Wu-Xu": [({}, [3])],
         "transition": [
             ({"edge_orders": [1, 1, 1], "variant": "equispaced"}, range(1, 6)),
@@ -98,7 +99,8 @@
         "Bernstein": [({}, range(1, 3))],
         "Taylor": [({}, range(0, 5))],
         "Bernardi-Raugel": [({}, [1, 2])],
-        "Guzman-Neilan": [({}, [1, 2])],
+        "Guzman-Neilan first kind": [({}, [1, 2])],
+        "Guzman-Neilan second kind": [({}, [1, 2])],
         "Wu-Xu": [({}, [4])],
         "transition": [
             (