Add DEFAULT_DEGREE for finite elements

FIAT/alfeld_sorokina.py

@@ -76,7 +76,10 @@ class AlfeldSorokina(finite_element.CiarletElement):
 
     This element belongs to a Stokes complex, and is paired with CG1(Alfeld).
     """
+    DEFAULT_DEGREE = 2
+
     def __init__(self, ref_el, degree=2):
+        degree = self._parse_degree(degree)
         dual = AlfeldSorokinaDualSet(ref_el, degree)
         poly_set = AlfeldSorokinaSpace(ref_el, degree)
         formdegree = ref_el.get_spatial_dimension() - 1  # (n-1)-form

FIAT/argyris.py

@@ -102,9 +102,10 @@ class Argyris(finite_element.CiarletElement):
     "integral(q)" -> dofs are evaluated by quadrature rules with the minimum
     degree required for unisolvence plus q.
     """
+    DEFAULT_DEGREE = 5
 
     def __init__(self, ref_el, degree=5, variant=None, quad_scheme=None):
-
+        degree = self._parse_degree(degree)
         splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             raise NotImplementedError(f"{type(self).__name__} is not implemented as a macroelement.")

FIAT/arnold_qin.py

@@ -58,7 +58,10 @@ def ArnoldQinSpace(ref_el, degree, reduced=False):
 class ArnoldQin(finite_element.CiarletElement):
     """The Arnold-Qin C0(Alfeld) quadratic macroelement with divergence in P0.
     This element belongs to a Stokes complex, and is paired with unsplit DG0."""
+    DEFAULT_DEGREE = 2
+
     def __init__(self, ref_el, degree=2, reduced=False):
+        degree = self._parse_degree(degree)
         poly_set = ArnoldQinSpace(ref_el, degree)
         if reduced:
             order = 1

FIAT/arnold_winther.py

@@ -123,7 +123,10 @@ def __init__(self, ref_el, degree=3):
 class ArnoldWinther(finite_element.CiarletElement):
     """The definition of the conforming Arnold-Winther element.
     """
+    DEFAULT_DEGREE = 3
+
     def __init__(self, ref_el, degree=3):
+        degree = self._parse_degree(degree)
         if ref_el.shape != TRIANGLE:
             raise ValueError(f"{type(self).__name__} only defined on triangles")
         Ps = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree)

FIAT/bell.py

@@ -48,10 +48,12 @@ def __init__(self, ref_el, degree):
 
 class Bell(finite_element.CiarletElement):
     """The Bell finite element."""
+    DEFAULT_DEGREE = 5
 
     def __init__(self, ref_el, degree=5):
         if ref_el.get_shape() != TRIANGLE:
             raise ValueError(f"{type(self).__name__} only defined on triangles")
+        degree = self._parse_degree(degree)
         if degree != 5:
             raise ValueError(f"{type(self).__name__} only defined for degree = 5.")
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)

FIAT/bernstein.py

@@ -50,8 +50,10 @@ def __init__(self, ref_el, degree):
 
 class Bernstein(FiniteElement):
     """A finite element with Bernstein polynomials as basis functions."""
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         dual = BernsteinDualSet(ref_el, degree)
         k = 0  # 0-form
         super().__init__(ref_el, dual, degree, k)

FIAT/brezzi_douglas_marini.py

@@ -90,16 +90,15 @@ class BrezziDouglasMarini(finite_element.CiarletElement):
     exactly. This is important when you want to have (nearly) div-preserving
     interpolation.
     """
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
 
         splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             ref_el = splitting(ref_el)
 
-        if degree < 1:
-            raise Exception("BDM_k elements only valid for k >= 1")
-
         sd = ref_el.get_spatial_dimension()
         if ref_el.is_macrocell():
             base_element = type(self)(ref_el.get_parent(), degree)

FIAT/c2_elements.py

@@ -85,7 +85,10 @@ def __init__(self, ref_complex, degree, vorder=None, reduced=False, quad_scheme=
 
 class BrambleZlamalC2(finite_element.CiarletElement):
     """The Bramble-Zlamal C2 element."""
+    DEFAULT_DEGREE = 9
+
     def __init__(self, ref_el, degree=9, reduced=False, quad_scheme=None):
+        degree = self._parse_degree(degree)
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
         dual = C2DualSet(ref_el, degree, reduced=reduced, quad_scheme=quad_scheme)
         super().__init__(poly_set, dual, degree, formdegree=0)
@@ -109,7 +112,10 @@ class AlfeldC2(finite_element.CiarletElement):
     """The Alfeld C^2 macroelement on a double barycentric split.
     See Section 7.5 of Lai & Schumacher for the quintic C^2 spline.
     """
+    DEFAULT_DEGREE = 5
+
     def __init__(self, ref_el, degree=5, reduced=False, quad_scheme=None):
+        degree = self._parse_degree(degree)
         poly_set = AlfeldC2Space(ref_el, degree)
         ref_complex = poly_set.get_reference_element()
         dual = C2DualSet(ref_complex, degree, reduced=reduced, quad_scheme=quad_scheme)

FIAT/crouzeix_raviart.py

@@ -76,8 +76,10 @@ class CrouzeixRaviart(finite_element.CiarletElement):
     Polynomial space:  P_k
     Dual basis:        Evaluation at points or integral moments
     """
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
         if degree % 2 != 1:
             raise ValueError("Crouzeix-Raviart only defined for odd degree")
 

FIAT/discontinuous_lagrange.py

@@ -214,6 +214,9 @@ class DiscontinuousLagrange(finite_element.CiarletElement):
                            variant='Alfeld' can be used to obtain a barycentrically refined
                               macroelement for Scott-Vogelius.
     """
+
+    DEFAULT_DEGREE = 0
+
     def __new__(cls, ref_el, degree, variant="equispaced"):
         if degree == 0:
             splitting, _ = parse_lagrange_variant(variant, discontinuous=True)
@@ -223,6 +226,8 @@ def __new__(cls, ref_el, degree, variant="equispaced"):
         return super().__new__(cls)
 
     def __init__(self, ref_el, degree, variant="equispaced"):
+        degree = self._parse_degree(degree)
+
         splitting, point_variant = parse_lagrange_variant(variant, discontinuous=True)
         if splitting is not None:
             ref_el = splitting(ref_el)

FIAT/finite_element.py

@@ -17,12 +17,18 @@
 from FIAT.quadrature_schemes import create_quadrature
 
 
-class FiniteElement(object):
+class FiniteElement:
     """Class implementing a basic abstraction template for general
     finite element families. Finite elements which inherit from
     this class are non-nodal unless they are CiarletElement subclasses.
     """
 
+    @property
+    def DEFAULT_DEGREE(self):
+        raise NotImplementedError(f"{type(self).__name__} does not specify a "
+                                  "default degree, please pass one explicitly "
+                                  "instead")
+
     def __init__(self, ref_el, dual, order, formdegree=None, mapping="affine", ref_complex=None):
         # Relevant attributes that do not necessarily depend on a PolynomialSet object:
         # The order (degree) of the polynomial basis
@@ -120,6 +126,12 @@ def is_nodal():
     def is_macroelement(self):
         return self.ref_el is not self.ref_complex
 
+    def _parse_degree(self, degree):
+        degree = self.DEFAULT_DEGREE if degree is None else degree
+        if degree < self.DEFAULT_DEGREE:
+            raise ValueError(f"{type(self).__name__} is only defined for degree >= {self.DEFAULT_DEGREE}.")
+        return degree
+
 
 class CiarletElement(FiniteElement):
     """Class implementing Ciarlet's abstraction of a finite element

FIAT/gopalakrishnan_lederer_schoberl.py

@@ -52,7 +52,10 @@ class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement):
     the weak symmetry constraint.
 
     """
+    DEFAULT_DEGREE = 0
+
     def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
 
         splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         assert variant == "integral"
@@ -81,6 +84,8 @@ def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree, variant=None, quad_sc
 
     Reference: https://doi.org/10.1093/imanum/drz022
     """
+    if degree < 1:
+        raise ValueError("GopalakrishnanLedererSchoberlFirstKind is only defined for degree >= 1.")
     fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree, variant=variant, quad_scheme=quad_scheme)
     entity_dofs = fe.entity_dofs()
     sd = ref_el.get_spatial_dimension()

FIAT/hct.py

@@ -80,7 +80,10 @@ class HsiehCloughTocher(finite_element.CiarletElement):
     super-smooth C^1 space from Groselj and Knez (2022) on a barycentric split,
     although there the basis functions are positive on an incenter split.
     """
+    DEFAULT_DEGREE = 3
+
     def __init__(self, ref_el, degree=3, reduced=False, quad_scheme=None):
+        degree = self._parse_degree(degree)
         ref_complex = macro.AlfeldSplit(ref_el)
         dual = HCTDualSet(ref_complex, degree, reduced=reduced, quad_scheme=quad_scheme)
         poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1, vorder=degree-1, variant="bubble")

FIAT/hellan_herrmann_johnson.py

@@ -99,9 +99,10 @@ class HellanHerrmannJohnson(finite_element.CiarletElement):
        HHJ(k) is the space of symmetric-matrix-valued polynomials of degree k
        or less with normal-normal continuity.
     """
-    def __init__(self, ref_el, degree=0, variant=None, quad_scheme=None):
-        if degree < 0:
-            raise ValueError(f"{type(self).__name__} only defined for degree >= 0")
+    DEFAULT_DEGREE = 0
+
+    def __init__(self, ref_el, degree=None, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
 
         splitting, variant, qdegree = check_format_variant(variant, degree)
         if splitting is not None:

FIAT/hermite.py

@@ -69,10 +69,12 @@ def __init__(self, ref_el):
 
 class CubicHermite(finite_element.CiarletElement):
     """The cubic Hermite finite element.  It is what it is."""
+    DEFAULT_DEGREE = 3
 
-    def __init__(self, ref_el, deg=3):
-        assert deg == 3
-        poly_set = polynomial_set.ONPolynomialSet(ref_el, 3)
+    def __init__(self, ref_el, degree=3):
+        degree = self._parse_degree(degree)
+        assert degree == 3
+        poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
         dual = CubicHermiteDualSet(ref_el)
 
-        super().__init__(poly_set, dual, 3)
+        super().__init__(poly_set, dual, degree)

FIAT/hierarchical.py

@@ -59,15 +59,18 @@ def __init__(self, ref_el, degree, codim=0, interpolant_deg=None, quad_scheme=No
 
 class Legendre(finite_element.CiarletElement):
     """Simplicial discontinuous element with Legendre polynomials."""
+    DEFAULT_DEGREE = 0
+
     def __new__(cls, ref_el, degree, variant=None):
-        if degree == 0:
+        if degree is None or degree == 0:
             splitting, variant, interpolant_deg = check_format_variant(variant, degree)
             if splitting is None and interpolant_deg == 0:
                 # FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size()
                 return P0(ref_el)
         return super().__new__(cls)
 
     def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
         splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             ref_el = splitting(ref_el)
@@ -102,12 +105,13 @@ def __init__(self, ref_el, degree, interpolant_deg=None, quad_scheme=None):
 
 class IntegratedLegendre(finite_element.CiarletElement):
     """Simplicial continuous element with integrated Legendre polynomials."""
+    DEFAULT_DEGREE = 1
+
     def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
         splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             ref_el = splitting(ref_el)
-        if degree < 1:
-            raise ValueError(f"{type(self).__name__} elements only valid for k >= 1")
         poly_set = ONPolynomialSet(ref_el, degree, variant="bubble")
         dual = IntegratedLegendreDual(ref_el, degree, interpolant_deg=interpolant_deg, quad_scheme=quad_scheme)
         formdegree = 0  # 0-form

FIAT/histopolation.py

@@ -60,7 +60,10 @@ def __init__(self, ref_el, degree):
 
 class Histopolation(finite_element.CiarletElement):
     """1D discontinuous element with integral DOFs on GLL subgrid."""
+    DEFAULT_DEGREE = 0
+
     def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         if ref_el.shape != LINE:
             raise ValueError("Histopolation elements are only defined in one dimension.")
 

FIAT/hu_zhang.py

@@ -85,9 +85,10 @@ def __init__(self, ref_el, degree, variant, qdegree, quad_scheme):
 
 class HuZhang(finite_element.CiarletElement):
     """The definition of the Hu-Zhang element."""
+    DEFAULT_DEGREE = 3
+
     def __init__(self, ref_el, degree=3, variant=None, quad_scheme=None):
-        if degree < 3:
-            raise ValueError(f"{type(self).__name__} only defined for degree >= 3")
+        degree = self._parse_degree(degree)
         if ref_el.shape != TRIANGLE:
             raise ValueError(f"{type(self).__name__} only defined on triangles")
         splitting, variant, qdegree = check_format_variant(variant, degree)

FIAT/johnson_mercier.py

@@ -58,8 +58,10 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None):
 
 class JohnsonMercier(finite_element.CiarletElement):
     """The Johnson-Mercier finite element."""
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
         ref_complex = macro.AlfeldSplit(ref_el)
         poly_set = macro.HDivSymPolynomialSet(ref_complex, degree)
         dual = JohnsonMercierDualSet(ref_complex, degree, variant=variant, quad_scheme=quad_scheme)

FIAT/kong_mulder_veldhuizen.py

@@ -86,8 +86,10 @@ class KongMulderVeldhuizen(finite_element.CiarletElement):
     W. A. MULDER
 
     """
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree, variant=None):
+        degree = self._parse_degree(degree)
         splitting, variant = parse_lagrange_variant(variant)
         if splitting:
             ref_el = splitting(ref_el)

FIAT/lagrange.py

@@ -72,7 +72,11 @@ class Lagrange(finite_element.CiarletElement):
                         entity id. The DOFs are always sorted by the entity ordering
                         and then lexicographically by lattice multiindex.
     """
+    DEFAULT_DEGREE = 1
+
     def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False):
+        degree = self._parse_degree(degree)
+
         splitting, point_variant = parse_lagrange_variant(variant)
         if splitting is not None:
             ref_el = splitting(ref_el)

FIAT/mixed.py

@@ -23,6 +23,9 @@ class MixedElement(FiniteElement):
 
     This object offers tabulation of the concatenated basis function
     tables along with an entity_dofs dict."""
+
+    DEFAULT_DEGREE = None
+
     def __init__(self, elements, ref_el=None):
         elements = tuple(elements)
 

FIAT/morley.py

@@ -53,10 +53,12 @@ def duals(ref_el, dim, degree):
 
 class Morley(finite_element.CiarletElement):
     """The Morley finite element."""
+    DEFAULT_DEGREE = 2
 
     def __init__(self, ref_el, degree=2):
         if ref_el.get_shape() not in {TRIANGLE, TETRAHEDRON}:
             raise ValueError("Morley only defined on simplices of dimension >= 2")
+        degree = self._parse_degree(degree)
         if degree != 2:
             raise ValueError("{type(self).__name__} only defined for degree == 2")
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)

FIAT/nedelec.py

@@ -197,8 +197,11 @@ class Nedelec(finite_element.CiarletElement):
     exactly. This is important when you want to have (nearly) curl-preserving
     interpolation.
     """
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
+
         splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             ref_el = splitting(ref_el)

FIAT/nedelec_second_kind.py

@@ -190,15 +190,15 @@ class NedelecSecondKind(CiarletElement):
     exactly. This is important when you want to have (nearly) curl-preserving
     interpolation.
     """
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
+        degree = self._parse_degree(degree)
+
         splitting, variant, interpolant_deg = check_format_variant(variant, degree)
         if splitting is not None:
             ref_el = splitting(ref_el)
 
-        # Check degree
-        assert degree >= 1, "Second kind Nedelecs start at 1!"
-
         # Get dimension
         d = ref_el.get_spatial_dimension()