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

@@ -30,10 +30,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:
@@ -61,7 +60,7 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False):
 
     default = "integral" if integral else "spectral"
     if integral:
-        supported_point_variants = {"integral": None, "fdm": "fdm", "demkowicz": "demkowicz"}
+        supported_point_variants = {"integral": None, "fdm": "fdm", "demkowicz": "demkowicz", "demkowicz-mass": "demkowicz-mass"}
     elif discontinuous:
         supported_point_variants = supported_dg_variants
     else:
@@ -71,7 +70,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:

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()
@@ -286,6 +293,7 @@ def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind):
     kind = 1
     variant = "fdm"
     variant = "demkowicz"
+    variant = "demkowicz-mass"
     # 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:

test/unit/test_hct.py

@@ -5,6 +5,8 @@
 from FIAT.reference_element import ufc_simplex
 from FIAT.functional import PointEvaluation
 from FIAT.macro import CkPolynomialSet
+from FIAT.quadrature_schemes import create_quadrature
+from FIAT.jacobi import eval_jacobi
 
 
 @pytest.fixture
@@ -74,3 +76,23 @@ def test_full_polynomials(cell, reduced):
     C1 = CkPolynomialSet(ref_complex, degree, order=1, variant="bubble")
     C1_tab = C1.tabulate(pts)[(0, 0)]
     assert span_greater_equal(tab, C1_tab)
+
+
+def test_reduced_normal_derivative(cell):
+    fe = HCT(cell, reduced=True)
+
+    ref_line = cell.construct_subelement(1)
+    Q = create_quadrature(ref_line, fe.degree()+2)
+    qpts, qwts = Q.get_points(), Q.get_weights()
+
+    bary = ref_line.compute_barycentric_coordinates(qpts)
+    leg2 = eval_jacobi(0, 0, 2, bary[:, 1] - bary[:, 0])
+    wts = numpy.multiply(leg2, qwts)
+    top = cell.get_topology()
+    for e in top[1]:
+        n = cell.compute_normal(e)
+        vals = fe.tabulate(1, qpts, entity=(1, e))
+        fn = vals[(1, 0)] * n[0] + vals[(0, 1)] * n[1]
+
+        fn = fn[:-3]
+        assert numpy.allclose(numpy.dot(fn, wts), 0)