Add BDDF element (#44)

* remove print

* correct docstring

* Add BDDM element on hex

Author: mscroggs

symfem/__init__.py

@@ -99,7 +99,7 @@ def create_element(cell_type, element_type, order):
         Raviart-Thomas, RT, N1div, dQ, NCE, RTCE, Qcurl, Q, NCF, RTCF, Qdiv, vector Q,
         vQ, Brezzi-Douglas-Marini, BDM, N2div, Morley, Hermite, Mardal-Tai-Winther,
         MTW, Argyris, bubble, dual, Buffa-Christiansen, BC, rotated Buffa-Christiansen,
-        RBC, Brezzi-Douglas-Fortin-Marini, BDFM
+        RBC, Brezzi-Douglas-Fortin-Marini, BDFM, Brezzi-Douglas-Duran-Fortin, BDDF
     order : int
         The order of the element.
     """

symfem/elements/bddm.py

@@ -0,0 +1,47 @@
+"""Brezzi-Douglas-Duran-Fortin elements."""
+
+from ..core.finite_element import FiniteElement
+from ..core.moments import make_integral_moment_dofs
+from ..core.polynomials import polynomial_set
+from ..core.symbolic import x, zero
+from ..core.calculus import curl
+from ..core.functionals import NormalIntegralMoment, IntegralMoment
+from .lagrange import DiscontinuousLagrange, VectorDiscontinuousLagrange
+
+
+def bddf_polyset(reference, order):
+    """Create the polynomial basis for a BDDF element."""
+    dim = reference.tdim
+    pset = []
+    assert reference.name == "hexahedron"
+    pset = polynomial_set(dim, dim, order)
+    pset.append(curl((zero, zero, x[0] ** (order + 1) * x[1])))
+    pset.append(curl((zero, x[0] * x[2] ** (order + 1), zero)))
+    pset.append(curl((x[1] ** (order + 1) * x[2], zero, zero)))
+    for i in range(1, order + 1):
+        pset.append(curl((zero, zero, x[0] * x[1] ** (i + 1) * x[2] ** (order - i))))
+        pset.append(curl((zero, x[0] ** (i + 1) * x[1] ** (order - i) * x[2], zero)))
+        pset.append(curl((x[0] ** (order - i) * x[1] * x[2] ** (i + 1), zero, zero)))
+
+    return pset
+
+
+class BDDF(FiniteElement):
+    """Brezzi-Douglas-Duran-Fortin Hdiv finite element."""
+
+    def __init__(self, reference, order):
+        poly = bddf_polyset(reference, order)
+
+        dofs = make_integral_moment_dofs(
+            reference,
+            facets=(NormalIntegralMoment, DiscontinuousLagrange, order),
+            cells=(IntegralMoment, VectorDiscontinuousLagrange, order - 2),
+        )
+
+        super().__init__(reference, order, poly, dofs, reference.tdim, reference.tdim)
+
+    names = ["Brezzi-Douglas-Duran-Fortin", "BDDF"]
+    references = ["hexahedron"]
+    min_order = 1
+    mapping = "contravariant"
+    continuity = "H(div)"

symfem/elements/bdfm.py

@@ -1,4 +1,4 @@
-"""Brezzi-Douglas-Fortin-Marini and Brezzi-Douglas-Duran-Fortin elements."""
+"""Brezzi-Douglas-Fortin-Marini elements."""
 
 from ..core.finite_element import FiniteElement
 from ..core.moments import make_integral_moment_dofs

test/utils.py

@@ -33,6 +33,7 @@
         "serendipity Hcurl": 3,
     },
     "hexahedron": {
+        "Brezzi-Douglas-Duran-Fortin": 2,
         "Brezzi-Douglas-Fortin-Marini": 2,
         "NCE": 2,
         "NCF": 2,
@@ -77,7 +78,6 @@ def elements(max_order=5, include_dual=True, include_non_dual=True,
                     max_o = max_o[r]
             else:
                 max_o = 100
-            print(r, e.names[0])
             if r in max_orders_getting_basis and e.names[0] in max_orders_getting_basis[r]:
                 max_o = min(max_orders_getting_basis[r][e.names[0]], max_o)