Simplex FDM variants for the De Rham complex

FIAT/brezzi_douglas_marini.py

@@ -6,7 +6,7 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 from FIAT import (finite_element, functional, dual_set,
-                  polynomial_set, nedelec)
+                  polynomial_set, nedelec, demkowicz)
 from FIAT.check_format_variant import check_format_variant
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.quadrature import FacetQuadratureRule
@@ -100,10 +100,15 @@ def __init__(self, ref_el, degree, variant=None):
             ref_el = splitting(ref_el)
 
         if degree < 1:
-            raise Exception("BDM_k elements only valid for k >= 1")
+            raise ValueError(f"{type(self).__name__} elements only valid for k >= 1")
 
         sd = ref_el.get_spatial_dimension()
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, ))
-        dual = BDMDualSet(ref_el, degree, variant, interpolant_deg)
+        if variant and variant.startswith("demkowicz"):
+            dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv", kind=2, variant=variant)
+        elif variant == "fdm":
+            dual = demkowicz.FDMDual(ref_el, degree, "HDiv", type(self))
+        else:
+            dual = BDMDualSet(ref_el, degree, variant, interpolant_deg)
         formdegree = sd - 1  # (n-1)-form
         super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")

FIAT/check_format_variant.py

@@ -32,10 +32,9 @@ def check_format_variant(variant, degree):
     if variant is None:
         variant = "integral"
 
-    match = re.match(r"^integral(?:\((\d+)\))?$", variant)
+    match = re.match(r"^(\w+)(?:\((\d+)\))?$", variant)
     if match:
-        variant = "integral"
-        extra_degree, = match.groups()
+        variant, extra_degree = match.groups()
         extra_degree = int(extra_degree) if extra_degree is not None else 0
         interpolant_degree = degree + extra_degree
         if interpolant_degree < degree:
@@ -63,7 +62,7 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False):
 
     default = "integral" if integral else "spectral"
     if integral:
-        supported_point_variants = {"integral": None, "point": "point"}
+        supported_point_variants = {"integral": None, "point": "point", "fdm": "fdm", "demkowicz": "demkowicz", "demkowiczmass": "demkowiczmass"}
     elif discontinuous:
         supported_point_variants = supported_dg_variants
     else:
@@ -73,7 +72,6 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False):
     splitting = None
     splitting_args = tuple()
     point_variant = supported_point_variants[default]
-
     for pre_opt in options:
         opt = pre_opt.lower()
         if opt in supported_splits: