mass decoupling variant

FIAT/brezzi_douglas_marini.py

@@ -94,8 +94,8 @@ def __init__(self, ref_el, degree, variant=None):
 
         sd = ref_el.get_spatial_dimension()
         poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, ))
-        if variant == "demkowicz":
-            dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv")
+        if variant and variant.startswith("demkowicz"):
+            dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv", variant=variant)
         elif variant == "fdm":
             dual = demkowicz.FDMDual(ref_el, degree, "HDiv", type(self))
         else:

FIAT/demkowicz.py

@@ -65,7 +65,7 @@ def map_duals(ref_el, dim, entity, mapping, Q_ref, Phis):
 
 class DemkowiczDual(DualSet):
 
-    def __init__(self, ref_el, degree, sobolev_space, kind=None):
+    def __init__(self, ref_el, degree, sobolev_space, kind=None, variant=None):
         nodes = []
         entity_ids = {}
         reduced_dofs = {}
@@ -76,6 +76,8 @@ def __init__(self, ref_el, degree, sobolev_space, kind=None):
         dual_mapping = {"HCurl": "contravariant", "HDiv": "covariant"}.get(sobolev_space, None)
         if kind is None:
             kind = 1 if formdegree == 0 else 2
+        if variant is None:
+            variant = "demkowicz"
 
         for dim in sorted(top):
             entity_ids[dim] = {}
@@ -91,7 +93,7 @@ def __init__(self, ref_el, degree, sobolev_space, kind=None):
                     entity_ids[dim][entity] = list(range(cur, len(nodes)))
                 reduced_dofs[dim] = len(nodes)
             else:
-                Q_ref, Phis, rdofs = self._reference_duals(dim, degree, formdegree, sobolev_space, kind)
+                Q_ref, Phis, rdofs = self._reference_duals(dim, degree, formdegree, sobolev_space, kind, variant)
                 reduced_dofs[dim] = rdofs
                 mapping = dual_mapping if dim == sd else trace
                 for entity in sorted(top[dim]):
@@ -101,9 +103,9 @@ def __init__(self, ref_el, degree, sobolev_space, kind=None):
                     entity_ids[dim][entity] = list(range(cur, len(nodes)))
 
         self._reduced_dofs = reduced_dofs
-        super(DemkowiczDual, self).__init__(nodes, ref_el, entity_ids)
+        super().__init__(nodes, ref_el, entity_ids)
 
-    def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind):
+    def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind, variant):
         facet = symmetric_simplex(dim)
         Q = create_quadrature(facet, 2 * degree)
         Qpts, Qwts = Q.get_points(), Q.get_weights()
@@ -122,10 +124,15 @@ def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind):
         if formdegree >= dim:
             K = inner(trial[:1], trial, Qwts)
         else:
-            dtrial = exterior_derivative(P_at_qpts)
-            if dim == 2 and formdegree == 1 and sobolev_space == "HDiv":
-                dtrial = dtrial[:, None, :]
-            K = self._bubble_derivative_moments(facet, degree, formdegree, kind, Qpts, Qwts, dtrial)
+            if variant == "demkowicz":
+                deriv = True
+                dtrial = exterior_derivative(P_at_qpts)
+                if dim == 2 and formdegree == 1 and sobolev_space == "HDiv":
+                    dtrial = dtrial[:, None, :]
+            else:
+                deriv = False
+                dtrial = trial
+            K = self._bubble_derivative_moments(facet, degree, formdegree, kind, Qpts, Qwts, dtrial, deriv=deriv)
         reduced_dofs = K.shape[0]
 
         # Evaluate type-II degrees of freedom on P
@@ -220,7 +227,7 @@ def __init__(self, ref_el, degree, sobolev_space, element):
         self.Q = Q
         self.V0 = phis[(0,) * sd]
         self.V1 = exterior_derivative(phis)
-        super(FDMDual, self).__init__(ref_el, degree, sobolev_space, kind=None)
+        super().__init__(ref_el, degree, sobolev_space, kind=None)
 
     def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind):
         entity_dofs = self.fe.entity_dofs()
@@ -285,7 +292,7 @@ def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind):
 
     kind = 1
     variant = "fdm"
-    variant = "demkowicz"
+    variant = "demkowicz(1)"
     # variant = None
     space_dict = {"H1": (CG, grad),
                   "HCurl": (N1Curl if kind == 1 else N2Curl, curl),

FIAT/hierarchical.py

@@ -78,7 +78,7 @@ def __init__(self, ref_el, degree, variant=None):
         if splitting is not None:
             ref_el = splitting(ref_el)
         poly_set = ONPolynomialSet(ref_el, degree)
-        if variant == "demkowicz":
+        if variant and variant.startswith("demkowicz"):
             dual = demkowicz.DemkowiczDual(ref_el, degree, "L2")
         elif variant == "fdm":
             dual = demkowicz.FDMDual(ref_el, degree, "L2", type(self))
@@ -135,7 +135,7 @@ def __init__(self, ref_el, degree, variant=None):
         if degree < 1:
             raise ValueError(f"{type(self).__name__} elements only valid for k >= 1")
         poly_set = ONPolynomialSet(ref_el, degree, variant="bubble")
-        if variant == "demkowicz":
+        if variant and variant.startswith("demkowicz"):
             dual = demkowicz.DemkowiczDual(ref_el, degree, "H1")
         elif variant == "fdm":
             dual = demkowicz.FDMDual(ref_el, degree, "H1", type(self))

FIAT/nedelec.py

@@ -196,8 +196,8 @@ class Nedelec(finite_element.CiarletElement):
 
     def __init__(self, ref_el, degree, variant=None):
 
-        if variant == "demkowicz":
-            dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl", kind=1)
+        if variant and variant.startswith("demkowicz"):
+            dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl", kind=1, variant=variant)
         elif variant == "fdm":
             dual = demkowicz.FDMDual(ref_el, degree, "HCurl", type(self))
         else:

FIAT/nedelec_second_kind.py

@@ -198,8 +198,8 @@ def __init__(self, ref_el, degree, variant=None):
 
         sd = ref_el.get_spatial_dimension()
         poly_set = ONPolynomialSet(ref_el, degree, (sd, ), variant="bubble")
-        if variant == "demkowicz":
-            dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl")
+        if variant and variant.startswith("demkowicz"):
+            dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl", variant=variant)
         elif variant == "fdm":
             dual = demkowicz.FDMDual(ref_el, degree, "HCurl", type(self))
         else:

FIAT/raviart_thomas.py

@@ -141,8 +141,8 @@ class RaviartThomas(finite_element.CiarletElement):
 
     def __init__(self, ref_el, degree, variant=None):
 
-        if variant == "demkowicz":
-            dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv", kind=1)
+        if variant and variant.startswith("demkowicz"):
+            dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv", kind=1, variant=variant)
         elif variant == "fdm":
             dual = demkowicz.FDMDual(ref_el, degree, "HDiv", type(self))
         else: