From d72b5e186e6aeb46cc3a5ea94f86d79a5a8d5545 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Tue, 24 Mar 2026 10:43:02 +0000 Subject: [PATCH 1/8] Add DEFAULT_DEGREE for finite elements --- FIAT/discontinuous_lagrange.py | 5 +++++ FIAT/finite_element.py | 13 ++++++++++++- FIAT/lagrange.py | 4 ++++ 3 files changed, 21 insertions(+), 1 deletion(-) diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py index fd08f8fa3..06896123c 100644 --- a/FIAT/discontinuous_lagrange.py +++ b/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) diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index f8d4447ba..e0f69ef8b 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -17,13 +17,21 @@ 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): + order = self._parse_degree(order) + # Relevant attributes that do not necessarily depend on a PolynomialSet object: # The order (degree) of the polynomial basis self.order = order @@ -120,6 +128,9 @@ def is_nodal(): def is_macroelement(self): return self.ref_el is not self.ref_complex + def _parse_degree(self, degree): + return self.DEFAULT_DEGREE if degree is None else degree + class CiarletElement(FiniteElement): """Class implementing Ciarlet's abstraction of a finite element diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index cb605ec89..b172d48a3 100644 --- a/FIAT/lagrange.py +++ b/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) From c39d926304a54fb3842a853c3f719267f2e7d0b8 Mon Sep 17 00:00:00 2001 From: Connor Ward Date: Tue, 24 Mar 2026 10:50:25 +0000 Subject: [PATCH 2/8] Add more defaults --- FIAT/mixed.py | 3 +++ FIAT/quadrature_element.py | 2 ++ 2 files changed, 5 insertions(+) diff --git a/FIAT/mixed.py b/FIAT/mixed.py index 590121845..5b729ae9f 100644 --- a/FIAT/mixed.py +++ b/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) diff --git a/FIAT/quadrature_element.py b/FIAT/quadrature_element.py index d61eff908..9d9675672 100644 --- a/FIAT/quadrature_element.py +++ b/FIAT/quadrature_element.py @@ -19,6 +19,8 @@ class QuadratureElement(FiniteElement): """A set of quadrature points pretending to be a finite element.""" + DEFAULT_DEGREE = None + def __init__(self, ref_el, points, weights=None): # Create entity dofs. entity_dofs = {dim: {entity: [] for entity in entities} From 22cd9165d072d35df134c16f5335a7b5463f0ff1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 24 Mar 2026 11:56:13 +0000 Subject: [PATCH 3/8] more elements --- FIAT/alfeld_sorokina.py | 3 +++ FIAT/argyris.py | 3 ++- FIAT/arnold_qin.py | 3 +++ FIAT/arnold_winther.py | 3 +++ FIAT/bell.py | 2 ++ FIAT/bernstein.py | 2 ++ FIAT/brezzi_douglas_marini.py | 5 ++--- FIAT/c2_elements.py | 6 ++++++ FIAT/crouzeix_raviart.py | 2 ++ FIAT/finite_element.py | 7 ++++--- FIAT/gopalakrishnan_lederer_schoberl.py | 5 +++++ FIAT/hct.py | 3 +++ FIAT/hellan_herrmann_johnson.py | 7 ++++--- FIAT/hermite.py | 10 ++++++---- FIAT/hierarchical.py | 10 +++++++--- FIAT/histopolation.py | 3 +++ FIAT/hu_zhang.py | 5 +++-- FIAT/johnson_mercier.py | 2 ++ FIAT/kong_mulder_veldhuizen.py | 2 ++ FIAT/morley.py | 2 ++ FIAT/nedelec.py | 3 +++ FIAT/nedelec_second_kind.py | 6 +++--- FIAT/powell_sabin.py | 6 ++++++ FIAT/raviart_thomas.py | 3 +++ FIAT/regge.py | 7 ++++--- FIAT/walkington.py | 2 ++ FIAT/wuxu.py | 6 ++++++ 27 files changed, 93 insertions(+), 25 deletions(-) diff --git a/FIAT/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py index e09f26d17..91192324e 100644 --- a/FIAT/alfeld_sorokina.py +++ b/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 diff --git a/FIAT/argyris.py b/FIAT/argyris.py index 29bef2d1d..4ab480a0f 100644 --- a/FIAT/argyris.py +++ b/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.") diff --git a/FIAT/arnold_qin.py b/FIAT/arnold_qin.py index 1c5794da8..72ff0fa5f 100644 --- a/FIAT/arnold_qin.py +++ b/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 diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index e7af81783..bfdbd567d 100644 --- a/FIAT/arnold_winther.py +++ b/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) diff --git a/FIAT/bell.py b/FIAT/bell.py index 8fccef82e..e9b28759b 100644 --- a/FIAT/bell.py +++ b/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) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index f960c64f4..07d291f58 100644 --- a/FIAT/bernstein.py +++ b/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) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 7a5f8906f..62aa886dd 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/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) diff --git a/FIAT/c2_elements.py b/FIAT/c2_elements.py index 66a6f2eda..cc6c6489a 100644 --- a/FIAT/c2_elements.py +++ b/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) diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index 1fd990cc9..ac5c61d3e 100644 --- a/FIAT/crouzeix_raviart.py +++ b/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") diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index e0f69ef8b..db2368da1 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -30,8 +30,6 @@ def DEFAULT_DEGREE(self): "instead") def __init__(self, ref_el, dual, order, formdegree=None, mapping="affine", ref_complex=None): - order = self._parse_degree(order) - # Relevant attributes that do not necessarily depend on a PolynomialSet object: # The order (degree) of the polynomial basis self.order = order @@ -129,7 +127,10 @@ def is_macroelement(self): return self.ref_el is not self.ref_complex def _parse_degree(self, degree): - return self.DEFAULT_DEGREE if degree is None else 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): diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py index aaaea94e6..4a0ce558f 100644 --- a/FIAT/gopalakrishnan_lederer_schoberl.py +++ b/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() diff --git a/FIAT/hct.py b/FIAT/hct.py index e93fbcdc5..f0421174c 100644 --- a/FIAT/hct.py +++ b/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") diff --git a/FIAT/hellan_herrmann_johnson.py b/FIAT/hellan_herrmann_johnson.py index 61e193028..e6ca76aef 100644 --- a/FIAT/hellan_herrmann_johnson.py +++ b/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: diff --git a/FIAT/hermite.py b/FIAT/hermite.py index 79c83ff97..bae03d52f 100644 --- a/FIAT/hermite.py +++ b/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) diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index 8428778de..dfa8e391f 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -59,8 +59,10 @@ 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() @@ -68,6 +70,7 @@ def __new__(cls, ref_el, degree, variant=None): 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 diff --git a/FIAT/histopolation.py b/FIAT/histopolation.py index ff5f06965..ad26576e0 100644 --- a/FIAT/histopolation.py +++ b/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.") diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index d1b0b121f..095427ac8 100644 --- a/FIAT/hu_zhang.py +++ b/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) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index af70ff0a9..e11c6ef2f 100644 --- a/FIAT/johnson_mercier.py +++ b/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) diff --git a/FIAT/kong_mulder_veldhuizen.py b/FIAT/kong_mulder_veldhuizen.py index c21c90596..071d6972f 100644 --- a/FIAT/kong_mulder_veldhuizen.py +++ b/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) diff --git a/FIAT/morley.py b/FIAT/morley.py index af09c52c3..55b4d1ccd 100644 --- a/FIAT/morley.py +++ b/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) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index bdd0b138a..bfd9c483f 100644 --- a/FIAT/nedelec.py +++ b/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) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index e4433feba..2b0352147 100644 --- a/FIAT/nedelec_second_kind.py +++ b/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() diff --git a/FIAT/powell_sabin.py b/FIAT/powell_sabin.py index e9b7c4f3f..e82c13575 100644 --- a/FIAT/powell_sabin.py +++ b/FIAT/powell_sabin.py @@ -44,7 +44,10 @@ class QuadraticPowellSabin6(finite_element.CiarletElement): """The PS6 macroelement is a C^1 quadratic macroelement defined on the 6-way Powell-Sabin split of a triangle. """ + DEFAULT_DEGREE = 2 + def __init__(self, ref_el, degree=2): + degree = self._parse_degree(degree) if degree != 2: raise ValueError("PS6 only defined for degree = 2") ref_complex = macro.PowellSabinSplit(ref_el) @@ -96,7 +99,10 @@ class QuadraticPowellSabin12(finite_element.CiarletElement): """The PS12 macroelement is a C^1 quadratic macroelement defined on the 12-way Powell-Sabin split of a triangle. """ + DEFAULT_DEGREE = 2 + def __init__(self, ref_el, degree=2): + degree = self._parse_degree(degree) if degree != 2: raise ValueError("PS12 only defined for degree = 2") ref_complex = macro.PowellSabin12Split(ref_el) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 60e090e51..45cc896b9 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -141,8 +141,11 @@ class RaviartThomas(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) diff --git a/FIAT/regge.py b/FIAT/regge.py index 8a8243028..a7ee0e567 100644 --- a/FIAT/regge.py +++ b/FIAT/regge.py @@ -59,9 +59,10 @@ class Regge(finite_element.CiarletElement): REG(k) is the space of symmetric-matrix-valued polynomials of degree k or less with tangential-tangential 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: diff --git a/FIAT/walkington.py b/FIAT/walkington.py index 9ed0b7f96..785ad050b 100644 --- a/FIAT/walkington.py +++ b/FIAT/walkington.py @@ -99,10 +99,12 @@ def __init__(self, ref_el, degree): class Walkington(finite_element.CiarletElement): """The Walkington C1 macroelement.""" + DEFAULT_DEGREE = 5 def __init__(self, ref_el, degree=5): if ref_el.get_shape() != TETRAHEDRON: raise ValueError(f"{type(self).__name__} only defined on tetrahedron") + degree = self._parse_degree(degree) if degree != 5: raise ValueError(f"{type(self).__name__} only defined for degree=5.") diff --git a/FIAT/wuxu.py b/FIAT/wuxu.py index e9102eca3..c586c8520 100644 --- a/FIAT/wuxu.py +++ b/FIAT/wuxu.py @@ -155,7 +155,10 @@ def __init__(self, ref_el, degree): class WuXuRobustH3NC(finite_element.CiarletElement): """The Wu-Xu robust H3 nonconforming finite element""" + DEFAULT_DEGREE = 7 + def __init__(self, ref_el, degree=7): + degree = self._parse_degree(degree) poly_set = WuXuH3NCSpace(ref_el, robust=True) assert degree == poly_set.degree dual = WuXuRobustH3NCDualSet(ref_el, degree) @@ -164,7 +167,10 @@ def __init__(self, ref_el, degree=7): class WuXuH3NC(finite_element.CiarletElement): """The Wu-Xu H3 nonconforming finite element""" + DEFAULT_DEGREE = 4 + def __init__(self, ref_el, degree=4): + degree = self._parse_degree(degree) poly_set = WuXuH3NCSpace(ref_el) assert degree == poly_set.degree dual = WuXuH3NCDualSet(ref_el, degree) From 862262493ecb03983809f78bcd320b6dd25cf687 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 24 Mar 2026 12:09:11 +0000 Subject: [PATCH 4/8] fixes --- FIAT/discontinuous_lagrange.py | 3 +-- FIAT/mixed.py | 3 --- FIAT/quadrature_element.py | 3 --- 3 files changed, 1 insertion(+), 8 deletions(-) diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py index 06896123c..1e046ae28 100644 --- a/FIAT/discontinuous_lagrange.py +++ b/FIAT/discontinuous_lagrange.py @@ -214,11 +214,10 @@ 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() diff --git a/FIAT/mixed.py b/FIAT/mixed.py index 5b729ae9f..590121845 100644 --- a/FIAT/mixed.py +++ b/FIAT/mixed.py @@ -23,9 +23,6 @@ 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) diff --git a/FIAT/quadrature_element.py b/FIAT/quadrature_element.py index 9d9675672..fd5dbb0d2 100644 --- a/FIAT/quadrature_element.py +++ b/FIAT/quadrature_element.py @@ -18,9 +18,6 @@ class QuadratureElement(FiniteElement): """A set of quadrature points pretending to be a finite element.""" - - DEFAULT_DEGREE = None - def __init__(self, ref_el, points, weights=None): # Create entity dofs. entity_dofs = {dim: {entity: [] for entity in entities} From fe2a76ec94d793b55779253590047741aeeb8f50 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 24 Mar 2026 12:13:01 +0000 Subject: [PATCH 5/8] serendipity --- FIAT/Sminus.py | 18 +++++++++--------- FIAT/serendipity.py | 3 ++- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 862199cb6..7613ca5b1 100644 --- a/FIAT/Sminus.py +++ b/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 = selt._parse_degree(degree) flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() self.fdim = dim @@ -406,10 +406,10 @@ 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: @@ -460,10 +460,10 @@ def __init__(self, ref_el, degree): 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: diff --git a/FIAT/serendipity.py b/FIAT/serendipity.py index 8ef7a6b81..fc11d2b61 100644 --- a/FIAT/serendipity.py +++ b/FIAT/serendipity.py @@ -49,6 +49,7 @@ def _replace_numbers_with_symbols(polynomials): class Serendipity(FiniteElement): + DEFAULT_DEGREE = 1 def __new__(cls, ref_el, degree): dim = ref_el.get_spatial_dimension() @@ -61,7 +62,7 @@ def __new__(cls, ref_el, degree): return self def __init__(self, ref_el, degree): - + degree = self._parse_degree(degree) flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() flat_topology = flat_el.get_topology() From 13fb50369f63bd430fd6420809bba8701a83adc5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 24 Mar 2026 12:17:21 +0000 Subject: [PATCH 6/8] Sminus --- FIAT/Sminus.py | 18 ++++++------------ FIAT/SminusCurl.py | 12 ++++++------ FIAT/SminusDiv.py | 12 ++++++------ FIAT/brezzi_douglas_marini_cube.py | 11 ++++------- 4 files changed, 22 insertions(+), 31 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 7613ca5b1..98a419e6f 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -48,7 +48,7 @@ class TrimmedSerendipity(FiniteElement): DEFAULT_DEGREE = 1 def __init__(self, ref_el, degree, mapping): - degree = selt._parse_degree(degree) + degree = self._parse_degree(degree) flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() self.fdim = dim @@ -412,9 +412,8 @@ 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,14 +447,9 @@ 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") @@ -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") diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 1f5736738..4c6eb284a 100644 --- a/FIAT/SminusCurl.py +++ b/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 = {} diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 26900357d..252807fc3 100644 --- a/FIAT/SminusDiv.py +++ b/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 = {} diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index ab42c5037..98bc7592e 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/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() From 092195654f6b094efec299278d36930cfdebc211 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 24 Mar 2026 12:27:38 +0000 Subject: [PATCH 7/8] Apply suggestions from code review Co-authored-by: Pablo Brubeck --- FIAT/lagrange.py | 1 - FIAT/nedelec.py | 1 - FIAT/nedelec_second_kind.py | 1 - FIAT/raviart_thomas.py | 1 - 4 files changed, 4 deletions(-) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index b172d48a3..160f8a029 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -76,7 +76,6 @@ class Lagrange(finite_element.CiarletElement): 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) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index bfd9c483f..ad77e2ff0 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -201,7 +201,6 @@ class Nedelec(finite_element.CiarletElement): 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) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 2b0352147..82d044a9a 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -194,7 +194,6 @@ class NedelecSecondKind(CiarletElement): 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) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 45cc896b9..94e6dfdc7 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -145,7 +145,6 @@ class RaviartThomas(finite_element.CiarletElement): 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) From 895839cd9fe76d159117583b419c1b3f0fb81b75 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 24 Mar 2026 13:57:10 +0000 Subject: [PATCH 8/8] more elements --- FIAT/arnold_winther.py | 3 +++ FIAT/christiansen_hu.py | 3 +++ FIAT/discontinuous_lagrange.py | 1 - FIAT/discontinuous_pc.py | 2 ++ FIAT/discontinuous_raviart_thomas.py | 3 ++- FIAT/discontinuous_taylor.py | 16 ++++++++-------- FIAT/gauss_radau.py | 3 +++ 7 files changed, 21 insertions(+), 10 deletions(-) diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index bfdbd567d..61f095930 100644 --- a/FIAT/arnold_winther.py +++ b/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) diff --git a/FIAT/christiansen_hu.py b/FIAT/christiansen_hu.py index 64fe89659..269fadef5 100644 --- a/FIAT/christiansen_hu.py +++ b/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) diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py index 1e046ae28..b8c6e370c 100644 --- a/FIAT/discontinuous_lagrange.py +++ b/FIAT/discontinuous_lagrange.py @@ -226,7 +226,6 @@ def __new__(cls, ref_el, degree, variant="equispaced"): 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) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index dbb566906..dd3c4c9cb 100644 --- a/FIAT/discontinuous_pc.py +++ b/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) diff --git a/FIAT/discontinuous_raviart_thomas.py b/FIAT/discontinuous_raviart_thomas.py index 8e3444e94..8abf8eb71 100644 --- a/FIAT/discontinuous_raviart_thomas.py +++ b/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") diff --git a/FIAT/discontinuous_taylor.py b/FIAT/discontinuous_taylor.py index 1438d1d3b..78fca46ca 100644 --- a/FIAT/discontinuous_taylor.py +++ b/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) diff --git a/FIAT/gauss_radau.py b/FIAT/gauss_radau.py index 784369b2b..01c3a14de 100644 --- a/FIAT/gauss_radau.py +++ b/FIAT/gauss_radau.py @@ -27,9 +27,12 @@ def __init__(self, ref_el, degree, right=True): class GaussRadau(finite_element.CiarletElement): """1D discontinuous element with nodes at the Gauss-Radau points.""" + DEFAULT_DEGREE = 0 + def __init__(self, ref_el, degree): if ref_el.shape != LINE: raise ValueError("Gauss-Radau elements are only defined in one dimension.") + degree = self._parse_degree(degree) poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) dual = GaussRadauDualSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form