make elements use polynomial subdegree as canonical degree (#302)

* make elements use polynomial subdegree as canonical degree

* fix tests

* defelement branch

* ruff

* BC degree

* remove assert degree == 1

* last_updated

* change

* remove print

* remove assert

* max order for test

* nedelec degree in demo

* only run coveralls on push to main (it takes >2 hours to run)

Author: mscroggs

.github/workflows/coveralls.yml

@@ -1,7 +1,6 @@
 name: 🥼 Coveralls
 
 on:
-  pull_request:
   push:
     branches:
       - main

.github/workflows/defelement.yml

@@ -27,7 +27,7 @@ jobs:
         with:
           path: ./defelement
           repository: DefElement/DefElement
-          ref: main
+          ref: mscroggs/remove-symfem-maps
       - name: Install requirements
         run: |
           cd defelement

CHANGELOG_SINCE_LAST_VERSION.md

@@ -2,3 +2,4 @@
 - Implement component-wise integrals of VectorFunction and MatrixFunction
 - Implement grad of VectorFunction
 - Add Gopalakrishnan-Lederer-Schoberl element
+- Use polynomial subdegree to index elements wherever possible

demo/nedelec.py

@@ -15,7 +15,7 @@
 from symfem.symbols import x
 from symfem.utils import allequal
 
-element = symfem.create_element("triangle", "Nedelec1", 4)
+element = symfem.create_element("triangle", "Nedelec1", 3)
 polys = element.get_polynomial_basis()
 
 # Check that the first 20 polynomials in the polynomial basis are

symfem/elements/abf.py

@@ -45,7 +45,7 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         dofs: ListOfFunctionals = make_integral_moment_dofs(
             reference,
             edges=(NormalIntegralMoment, Lagrange, order, {"variant": variant}),
-            faces=(IntegralMoment, Nedelec, order, {"variant": variant}),
+            faces=(IntegralMoment, Nedelec, order - 1, {"variant": variant}),
         )
 
         for i in range(order + 1):
@@ -98,4 +98,4 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
     min_order = 0
     continuity = "H(div)"
     value_type = "vector"
-    last_updated = "2023.06"
+    last_updated = "2025.03"

symfem/elements/aw.py

@@ -40,25 +40,24 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         self.variant = variant
         poly: typing.List[FunctionInput] = []
         poly += [
-            ((p[0], p[1]), (p[1], p[2]))
-            for p in polynomial_set_vector(reference.tdim, 3, order - 1)
+            ((p[0], p[1]), (p[1], p[2])) for p in polynomial_set_vector(reference.tdim, 3, order)
         ]
         poly += [
             (
                 (
-                    (order - k + 1) * (order - k + 2) * x[0] ** k * x[1] ** (order - k),
-                    -k * (order - k + 2) * x[0] ** (k - 1) * x[1] ** (order - k + 1),
+                    (order - k + 2) * (order - k + 3) * x[0] ** k * x[1] ** (order - k + 1),
+                    -k * (order - k + 3) * x[0] ** (k - 1) * x[1] ** (order - k + 2),
                 ),
                 (
-                    -k * (order - k + 2) * x[0] ** (k - 1) * x[1] ** (order - k + 1),
-                    k * (k - 1) * x[0] ** (k - 2) * x[1] ** (order - k + 2),
+                    -k * (order - k + 3) * x[0] ** (k - 1) * x[1] ** (order - k + 2),
+                    k * (k - 1) * x[0] ** (k - 2) * x[1] ** (order - k + 3),
                 ),
             )
-            for k in range(order + 1)
+            for k in range(order + 2)
         ]
         poly += [
-            ((0, x[0] ** order), (x[0] ** order, -order * x[0] ** (order - 1) * x[1])),
-            ((0, 0), (0, x[0] ** order)),
+            ((0, x[0] ** (order + 1)), (x[0] ** (order + 1), -(order + 1) * x[0] ** order * x[1])),
+            ((0, 0), (0, x[0] ** (order + 1))),
         ]
 
         dofs: ListOfFunctionals = []
@@ -71,7 +70,7 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
                 )
         for e_n in range(reference.sub_entity_count(1)):
             sub_ref = reference.sub_entity(1, e_n)
-            sub_e = Lagrange(sub_ref.default_reference(), order - 2, variant)
+            sub_e = Lagrange(sub_ref.default_reference(), order - 1, variant)
             for dof_n, dof in enumerate(sub_e.dofs):
                 p = sub_e.get_basis_function(dof_n).get_function()
                 for component in [sub_ref.normal(), sub_ref.tangent()]:
@@ -86,7 +85,7 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
                             mapping="double_contravariant",
                         )
                     )
-        sub_e = Lagrange(reference, order - 3, variant)
+        sub_e = Lagrange(reference, order - 2, variant)
         for dof_n, dof in enumerate(sub_e.dofs):
             p = sub_e.get_basis_function(dof_n).get_function()
             for component22 in [((1, 0), (0, 0)), ((0, 1), (0, 0)), ((0, 0), (0, 1))]:
@@ -100,10 +99,10 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
                     )
                 )
 
-        if order >= 4:
-            sub_e = Lagrange(reference, order - 4, variant)
+        if order >= 3:
+            sub_e = Lagrange(reference, order - 3, variant)
             for p, dof in zip(sub_e.get_basis_functions(), sub_e.dofs):
-                if sympy.Poly(p.as_sympy(), x[:2]).degree() != order - 4:
+                if sympy.Poly(p.as_sympy(), x[:2]).degree() != order - 3:
                     continue
                 f = p * x[0] ** 2 * x[1] ** 2 * (1 - x[0] - x[1]) ** 2
                 J = [
@@ -134,26 +133,26 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
 
     @property
     def lagrange_subdegree(self) -> int:
-        return self.order - 1
+        return self.order
 
     @property
     def lagrange_superdegree(self) -> typing.Optional[int]:
-        return self.order
+        return self.order + 1
 
     @property
     def polynomial_subdegree(self) -> int:
-        return self.order - 1
+        return self.order
 
     @property
     def polynomial_superdegree(self) -> typing.Optional[int]:
-        return self.order
+        return self.order + 1
 
     names = ["Arnold-Winther", "AW", "conforming Arnold-Winther"]
     references = ["triangle"]
-    min_order = 3
+    min_order = 2
     continuity = "integral inner H(div)"
     value_type = "symmetric matrix"
-    last_updated = "2024.09"
+    last_updated = "2025.03"
 
 
 class NonConformingArnoldWinther(CiarletElement):
@@ -171,8 +170,7 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         self.variant = variant
         poly: typing.List[FunctionInput] = []
         poly += [
-            ((p[0], p[1]), (p[1], p[2]))
-            for p in polynomial_set_vector(reference.tdim, 3, order - 1)
+            ((p[0], p[1]), (p[1], p[2])) for p in polynomial_set_vector(reference.tdim, 3, order)
         ]
 
         poly += [
@@ -235,24 +233,24 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
 
     @property
     def lagrange_subdegree(self) -> int:
-        return self.order - 1
+        return self.order
 
     @property
     def lagrange_superdegree(self) -> typing.Optional[int]:
-        return self.order
+        return self.order + 1
 
     @property
     def polynomial_subdegree(self) -> int:
-        return self.order - 1
+        return self.order
 
     @property
     def polynomial_superdegree(self) -> typing.Optional[int]:
-        return self.order
+        return self.order + 1
 
     names = ["nonconforming Arnold-Winther", "nonconforming AW"]
     references = ["triangle"]
-    min_order = 2
-    max_order = 2
+    min_order = 1
+    max_order = 1
     continuity = "integral inner H(div)"
     value_type = "symmetric matrix"
-    last_updated = "2024.10"
+    last_updated = "2025.03"

symfem/elements/bdfm.py

@@ -33,38 +33,38 @@ def bdfm_polyset(reference: Reference, order: int) -> typing.List[FunctionInput]
     """
     dim = reference.tdim
     pset: typing.List[FunctionInput] = []
-    pset += polynomial_set_vector(dim, dim, order - 1)
+    pset += polynomial_set_vector(dim, dim, order)
     if reference.name == "quadrilateral":
-        for i in range(1, order + 1):
-            j = order - i
+        for i in range(1, order + 2):
+            j = order + 1 - i
             pset.append((x[0] ** i * x[1] ** j, 0))
             pset.append((0, x[1] ** i * x[0] ** j))
     elif reference.name == "triangle":
-        for i in range(order - 1):
-            p = x[0] ** i * x[1] ** (order - 2 - i)
+        for i in range(order):
+            p = x[0] ** i * x[1] ** (order - 1 - i)
             pset.append((x[0] * (x[0] + x[1]) * p, 0))
             pset.append((0, x[1] * (x[0] + x[1]) * p))
-        p = x[0] ** (order - 1)
+        p = x[0] ** order
         pset.append((x[0] * p, x[1] * p))
     elif reference.name == "hexahedron":
-        for i in range(1, order + 1):
-            for j in range(order + 1 - i):
-                k = order - i - j
+        for i in range(1, order + 2):
+            for j in range(order + 2 - i):
+                k = order + 1 - i - j
                 pset.append((x[0] ** i * x[1] ** j * x[2] ** k, 0, 0))
                 pset.append((0, x[1] ** i * x[0] ** j * x[2] ** k, 0))
                 pset.append((0, 0, x[2] ** i * x[0] ** j * x[1] ** k))
     elif reference.name == "tetrahedron":
-        for i in range(order - 1):
-            for j in range(order - i - 1):
-                p = x[0] ** i * x[1] ** j * x[2] ** (order - 2 - i - j)
+        for i in range(order):
+            for j in range(order - i):
+                p = x[0] ** i * x[1] ** j * x[2] ** (order - 1 - i - j)
                 pset.append((x[0] * (x[0] + x[1] + x[2]) * p, 0, 0))
                 pset.append((0, x[1] * (x[0] + x[1] + x[2]) * p, 0))
                 pset.append((0, 0, x[2] * (x[0] + x[1] + x[2]) * p))
-        for i in range(order):
-            p = x[0] ** i * x[1] ** (order - 1 - i)
+        for i in range(order + 1):
+            p = x[0] ** i * x[1] ** (order - i)
             pset.append((x[0] * p, x[1] * p, x[2] * p))
-        for i in range(1, order):
-            p = x[0] ** i * x[1] ** (order - 1 - i)
+        for i in range(1, order + 1):
+            p = x[0] ** i * x[1] ** (order - i)
             pset.append((p * (x[0] + x[1]), 0, x[1] * p))
 
     return pset
@@ -90,14 +90,14 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         if reference.name in ["triangle", "tetrahedron"]:
             dofs = make_integral_moment_dofs(
                 reference,
-                facets=(NormalIntegralMoment, Lagrange, order - 1, {"variant": variant}),
+                facets=(NormalIntegralMoment, Lagrange, order, {"variant": variant}),
                 cells=(IntegralMoment, NedelecFirstKind, order - 1, {"variant": variant}),
             )
         else:
             dofs = make_integral_moment_dofs(
                 reference,
-                facets=(NormalIntegralMoment, DPC, order - 1, {"variant": variant}),
-                cells=(IntegralMoment, VectorDPC, order - 2, {"variant": variant}),
+                facets=(NormalIntegralMoment, DPC, order, {"variant": variant}),
+                cells=(IntegralMoment, VectorDPC, order - 1, {"variant": variant}),
             )
         self.variant = variant
 
@@ -114,25 +114,25 @@ def init_kwargs(self) -> typing.Dict[str, typing.Any]:
     @property
     def lagrange_subdegree(self) -> int:
         if self.reference.name in ["triangle", "tetrahedron"]:
-            return self.order - 1
+            return self.order
         else:
-            return self.order // self.reference.tdim
+            return (self.order + 1) // self.reference.tdim
 
     @property
     def lagrange_superdegree(self) -> typing.Optional[int]:
-        return self.order
+        return self.order + 1
 
     @property
     def polynomial_subdegree(self) -> int:
-        return self.order - 1
+        return self.order
 
     @property
     def polynomial_superdegree(self) -> typing.Optional[int]:
-        return self.order
+        return self.order + 1
 
     names = ["Brezzi-Douglas-Fortin-Marini", "BDFM"]
     references = ["triangle", "quadrilateral", "hexahedron", "tetrahedron"]
-    min_order = 1
+    min_order = 0
     continuity = "H(div)"
     value_type = "vector"
-    last_updated = "2024.10"
+    last_updated = "2025.03"

symfem/elements/bdm.py

@@ -38,7 +38,7 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
         dofs: ListOfFunctionals = make_integral_moment_dofs(
             reference,
             facets=(NormalIntegralMoment, Lagrange, order, {"variant": variant}),
-            cells=(IntegralMoment, NedelecFirstKind, order - 1, {"variant": variant}),
+            cells=(IntegralMoment, NedelecFirstKind, order - 2, {"variant": variant}),
         )
         self.variant = variant
 
@@ -73,4 +73,4 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
     min_order = 1
     continuity = "H(div)"
     value_type = "vector"
-    last_updated = "2023.06"
+    last_updated = "2025.03"

symfem/elements/bell.py

@@ -27,7 +27,6 @@ def __init__(self, reference: Reference, order: int, variant: str = "equispaced"
             variant: The variant of the element
         """
         assert reference.name == "triangle"
-        assert order == 5
         dofs: ListOfFunctionals = []
         for v_n, v in enumerate(reference.vertices):
             dofs.append(PointEvaluation(reference, v, entity=(0, v_n)))
@@ -72,8 +71,9 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
 
     names = ["Bell"]
     references = ["triangle"]
-    min_order = 5
-    max_order = 5
+    min_order = 4
+    max_order = 4
     continuity = "C1"
     value_type = "scalar"
-    last_updated = "2023.05"
+    last_updated = "2025.03"
+    _max_continuity_test_order = 3

symfem/elements/dual.py

@@ -68,11 +68,6 @@ def get_polynomial_basis(self, reshape: bool = True) -> typing.List[Function]:
         """
         raise ValueError("Polynomial basis not supported for barycentric dual elements.")
 
-    @property
-    def maximum_degree(self) -> int:
-        """Get the maximum degree of this polynomial set for the element."""
-        raise NotImplementedError()
-
     def get_dual_matrix(self) -> sympy.matrices.dense.MutableDenseMatrix:
         """Get the dual matrix.
 
@@ -287,7 +282,6 @@ def __init__(self, reference: DualPolygon, order: int):
             reference: The reference element
             order: The polynomial order
         """
-        assert order == 1
         dual_coefficients: typing.List[
             typing.List[typing.List[typing.Union[int, sympy.core.expr.Expr]]]
         ] = [
@@ -339,11 +333,11 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
 
     names = ["Buffa-Christiansen", "BC"]
     references = ["dual polygon"]
-    min_order = 1
-    max_order = 1
+    min_order = 0
+    max_order = 0
     continuity = "H(div)"
     value_type = "vector dual"
-    last_updated = "2023.05"
+    last_updated = "2025.03"
 
 
 class RotatedBuffaChristiansen(DualCiarletElement):
@@ -356,7 +350,6 @@ def __init__(self, reference: DualPolygon, order: int):
             reference: The reference element
             order: The polynomial order
         """
-        assert order == 1
         dual_coefficients: typing.List[
             typing.List[typing.List[typing.Union[int, sympy.core.expr.Expr]]]
         ] = [
@@ -408,8 +401,8 @@ def polynomial_superdegree(self) -> typing.Optional[int]:
 
     names = ["rotated Buffa-Christiansen", "RBC"]
     references = ["dual polygon"]
-    min_order = 1
-    max_order = 1
+    min_order = 0
+    max_order = 0
     continuity = "H(div)"
     value_type = "vector dual"
-    last_updated = "2023.05"
+    last_updated = "2025.03"