Add DEFAULT_DEGREE for finite elements

FIAT/Sminus.py

@@ -45,10 +45,10 @@ def choose_ijk_total(degree):
 
 
 class TrimmedSerendipity(FiniteElement):
-    def __init__(self, ref_el, degree, mapping):
-        if degree < 1:
-            raise Exception("Trimmed serendipity elements only valid for k >= 1")
+    DEFAULT_DEGREE = 1
 
+    def __init__(self, ref_el, degree, mapping):
+        degree = self._parse_degree(degree)
         flat_el = flatten_reference_cube(ref_el)
         dim = flat_el.get_spatial_dimension()
         self.fdim = dim
@@ -406,15 +406,14 @@ def I_lambda_1_tilde_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid):
 
 # This is always 1-forms regardless of 2 or 3 dimensions.
 class TrimmedSerendipityEdge(TrimmedSerendipity):
-    def __init__(self, ref_el, degree):
-        if degree < 1:
-            raise Exception("Trimmed Serendipity_k edge elements only valid for k >= 1")
+    DEFAULT_DEGREE = 1
 
+    def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         flat_el = flatten_reference_cube(ref_el)
         dim = flat_el.get_spatial_dimension()
-        if dim != 2:
-            if dim != 3:
-                raise Exception("Trimmed Serendipity_k edge elements only valid for dimensions 2 and 3")
+        if dim not in {2, 3}:
+            raise Exception("Trimmed Serendipity_k edge elements only valid for dimensions 2 and 3")
 
         verts = flat_el.get_vertices()
 
@@ -448,22 +447,17 @@ def __init__(self, ref_el, degree):
             else:
                 IL = ()
 
-        Sminus_list = EL + FL
-        if dim == 3:
-            Sminus_list = Sminus_list + IL
+        Sminus_list = EL + FL + IL
 
-        if dim == 2:
-            self.basis = {(0, 0): Array(Sminus_list)}
-        else:
-            self.basis = {(0, 0, 0): Array(Sminus_list)}
+        self.basis = {(0,)*dim: Array(Sminus_list)}
         super().__init__(ref_el=ref_el, degree=degree, mapping="covariant piola")
 
 
 class TrimmedSerendipityFace(TrimmedSerendipity):
-    def __init__(self, ref_el, degree):
-        if degree < 1:
-            raise Exception("Trimmed serendipity face elements only valid for k >= 1")
+    DEFAULT_DEGREE = 1
 
+    def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         flat_el = flatten_reference_cube(ref_el)
         dim = flat_el.get_spatial_dimension()
         if dim != 2:
@@ -484,5 +478,5 @@ def __init__(self, ref_el, degree):
 
         Sminus_list = EL + FL
         Sminus_list = [[-a[1], a[0]] for a in Sminus_list]
-        self.basis = {(0, 0): Array(Sminus_list)}
+        self.basis = {(0,)*dim: Array(Sminus_list)}
         super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola")

FIAT/SminusCurl.py

@@ -26,16 +26,16 @@ def choose_ijk_total(degree):
 
 
 class TrimmedSerendipityCurl(FiniteElement):
-    def __init__(self, ref_el, degree):
-        if degree < 1:
-            raise Exception("Trimmed serendipity elements only valid for k >= 1")
 
+    DEFAULT_DEGREE = 1
+
+    def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         flat_el = flatten_reference_cube(ref_el)
         dim = flat_el.get_spatial_dimension()
         self.fdim = dim
-        if dim != 3:
-            if dim != 2:
-                raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3")
+        if dim not in {2, 3}:
+            raise ValueError(f"{type(self).__name__} only valid for dimensions 2 and 3")
 
         flat_topology = flat_el.get_topology()
         entity_ids = {}

FIAT/SminusDiv.py

@@ -26,16 +26,16 @@ def choose_ijk_total(degree):
 
 
 class TrimmedSerendipityDiv(FiniteElement):
-    def __init__(self, ref_el, degree):
-        if degree < 1:
-            raise Exception("Trimmed serendipity elements only valid for k >= 1")
 
+    DEFAULT_DEGREE = 1
+
+    def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         flat_el = flatten_reference_cube(ref_el)
         dim = flat_el.get_spatial_dimension()
         self.fdim = dim
-        if dim != 3:
-            if dim != 2:
-                raise Exception("Trimmed serendipity elements only valid for dimensions 2 and 3")
+        if dim not in {2, 3}:
+            raise ValueError(f"{type(self).__name__} only valid for dimensions 2 and 3")
 
         flat_topology = flat_el.get_topology()
         entity_ids = {}

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

@@ -62,7 +62,10 @@ def __init__(self, ref_el, degree=2):
 class ArnoldWintherNC(finite_element.CiarletElement):
     """The definition of the nonconforming Arnold-Winther element.
     """
+    DEFAULT_DEGREE = 2
+
     def __init__(self, ref_el, degree=2):
+        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)
@@ -123,7 +126,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/brezzi_douglas_marini_cube.py

@@ -39,22 +39,19 @@ class BrezziDouglasMariniCube(FiniteElement):
     The Brezzi-Douglas-Marini element on quadrilateral cells.
 
     :arg ref_el: The reference element.
-    :arg k: The degree.
+    :arg degree: The degree.
     :arg mapping: A string giving the Piola mapping.
                   Either 'contravariant Piola' or 'covariant Piola'.
     """
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree, mapping):
-
-        # Check that ref_el and degree are appropriate
-        if degree < 1:
-            raise Exception("BDMc_k elements only valid for k >= 1")
-
+        degree = self._parse_degree(degree)
         flat_el = flatten_reference_cube(ref_el)
         dim = flat_el.get_spatial_dimension()
         self.fdim = dim
         if dim != 2:
-            raise Exception("BDMc_k elements only valid for dimension 2")
+            raise Exception(f"{type(self).__name__} only valid for dimension 2")
 
         # Collect the IDs of the reference element entities
         flat_topology = flat_el.get_topology()

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/christiansen_hu.py

@@ -66,7 +66,10 @@ def ChristiansenHuSpace(ref_el, degree, reduced=False):
 class ChristiansenHu(finite_element.CiarletElement):
     """The Christiansen-Hu C^0(Worsey-Farin) linear macroelement with divergence in P0.
     This element belongs to a Stokes complex, and is paired with unsplit DG0."""
+    DEFAULT_DEGREE = 1
+
     def __init__(self, ref_el, degree=1):
+        degree = self._parse_degree(degree)
         if degree != 1:
             raise ValueError("Christiansen-Hu only defined for degree = 1")
         poly_set = ChristiansenHuSpace(ref_el, degree)

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,15 +214,18 @@ 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:
+        if degree is None or degree == 0:
             splitting, _ = parse_lagrange_variant(variant, discontinuous=True)
             if splitting is None and not ref_el.is_macrocell():
                 # FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size()
                 return P0.P0(ref_el)
         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/discontinuous_pc.py

@@ -99,8 +99,10 @@ def __init__(self, ref_el, flat_el, degree):
 
 class HigherOrderDPC(finite_element.CiarletElement):
     """The DPC finite element.  It is what it is."""
+    DEFAULT_DEGREE = 0
 
     def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         flat_el = flatten_reference_cube(ref_el)
         poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], degree)
         dual = DPCDualSet(ref_el, flat_el, degree)

FIAT/discontinuous_raviart_thomas.py

@@ -54,9 +54,10 @@ def __init__(self, ref_el, degree):
 
 class DiscontinuousRaviartThomas(finite_element.CiarletElement):
     """The discontinuous Raviart-Thomas finite element"""
+    DEFAULT_DEGREE = 1
 
     def __init__(self, ref_el, degree):
-
+        degree = self._parse_degree(degree)
         poly_set = RTSpace(ref_el, degree)
         dual = DRTDualSet(ref_el, degree)
         super().__init__(poly_set, dual, degree, mapping="contravariant piola")

FIAT/discontinuous_taylor.py

@@ -39,18 +39,18 @@ def __init__(self, ref_el, degree):
         super().__init__(nodes, ref_el, entity_ids)
 
 
-class HigherOrderDiscontinuousTaylor(finite_element.CiarletElement):
+class DiscontinuousTaylor(finite_element.CiarletElement):
     """The discontinuous Taylor finite element. Use a Taylor basis for DG."""
+    DEFAULT_DEGREE = 0
+
+    def __new__(cls, ref_el, degree):
+        if degree is None or degree == 0:
+            return P0.P0(ref_el)
+        return super().__new__(cls)
 
     def __init__(self, ref_el, degree):
+        degree = self._parse_degree(degree)
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)
         dual = DiscontinuousTaylorDualSet(ref_el, degree)
         formdegree = ref_el.get_spatial_dimension()  # n-form
         super().__init__(poly_set, dual, degree, formdegree)
-
-
-def DiscontinuousTaylor(ref_el, degree):
-    if degree == 0:
-        return P0.P0(ref_el)
-    else:
-        return HigherOrderDiscontinuousTaylor(ref_el, degree)

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