diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 862199cb..98a419e6 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 = 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") diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index 1f573673..4c6eb284 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 26900357..252807fc 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/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py index e09f26d1..91192324 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 29bef2d1..4ab480a0 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 1c5794da..72ff0fa5 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 e7af8178..61f09593 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) @@ -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) diff --git a/FIAT/bell.py b/FIAT/bell.py index 8fccef82..e9b28759 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 f960c64f..07d291f5 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 7a5f8906..62aa886d 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/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index ab42c503..98bc7592 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() diff --git a/FIAT/c2_elements.py b/FIAT/c2_elements.py index 66a6f2ed..cc6c6489 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/christiansen_hu.py b/FIAT/christiansen_hu.py index 64fe8965..269fadef 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/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index 1fd990cc..ac5c61d3 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/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py index fd08f8fa..b8c6e370 100644 --- a/FIAT/discontinuous_lagrange.py +++ b/FIAT/discontinuous_lagrange.py @@ -214,8 +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() @@ -223,6 +225,7 @@ 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/discontinuous_pc.py b/FIAT/discontinuous_pc.py index dbb56690..dd3c4c9c 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 8e3444e9..8abf8eb7 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 1438d1d3..78fca46c 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/finite_element.py b/FIAT/finite_element.py index f8d4447b..db2368da 100644 --- a/FIAT/finite_element.py +++ b/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 diff --git a/FIAT/gauss_radau.py b/FIAT/gauss_radau.py index 784369b2..01c3a14d 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 diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py index aaaea94e..4a0ce558 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 e93fbcdc..f0421174 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 61e19302..e6ca76ae 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 79c83ff9..bae03d52 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 8428778d..dfa8e391 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 ff5f0696..ad26576e 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 d1b0b121..095427ac 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 af70ff0a..e11c6ef2 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 c21c9059..071d6972 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/lagrange.py b/FIAT/lagrange.py index cb605ec8..160f8a02 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -72,7 +72,10 @@ 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) diff --git a/FIAT/morley.py b/FIAT/morley.py index af09c52c..55b4d1cc 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 bdd0b138..ad77e2ff 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -197,8 +197,10 @@ 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 e4433feb..82d044a9 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -190,15 +190,14 @@ 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 e9b7c4f3..e82c1357 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/quadrature_element.py b/FIAT/quadrature_element.py index d61eff90..fd5dbb0d 100644 --- a/FIAT/quadrature_element.py +++ b/FIAT/quadrature_element.py @@ -18,7 +18,6 @@ class QuadratureElement(FiniteElement): """A set of quadrature points pretending to be a finite element.""" - def __init__(self, ref_el, points, weights=None): # Create entity dofs. entity_dofs = {dim: {entity: [] for entity in entities} diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index 60e090e5..94e6dfdc 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -141,8 +141,10 @@ 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 8a824302..a7ee0e56 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/serendipity.py b/FIAT/serendipity.py index 8ef7a6b8..fc11d2b6 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() diff --git a/FIAT/walkington.py b/FIAT/walkington.py index 9ed0b7f9..785ad050 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 e9102eca..c586c852 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)