Merged in l2pullbacks (pull request #114)

L2pullbacks

Approved-by: David Ham <david.ham@imperial.ac.uk>

test/test_apply_function_pullbacks.py

@@ -34,6 +34,7 @@ def test_apply_single_function_pullbacks_triangle3d():
     cell = triangle3d
     domain = as_domain(cell)
 
+    UL2 = FiniteElement("DG L2", cell, 1)
     U0 = FiniteElement("DG", cell, 0)
     U = FiniteElement("CG", cell, 1)
     V = VectorElement("CG", cell, 1)
@@ -44,6 +45,7 @@ def test_apply_single_function_pullbacks_triangle3d():
     COV2T = FiniteElement("Regge", cell, 0)   # (0, 2)-symmetric tensors
     CONTRA2T = FiniteElement("HHJ", cell, 0)  # (2, 0)-symmetric tensors
 
+    Uml2 = UL2*UL2
     Um = U*U
     Vm = U*V
     Vdm = V*Vd
@@ -55,6 +57,7 @@ def test_apply_single_function_pullbacks_triangle3d():
 
     W = S*T*Vc*Vd*V*U
 
+    ul2 = Coefficient(UL2)
     u = Coefficient(U)
     v = Coefficient(V)
     vd = Coefficient(Vd)
@@ -64,6 +67,7 @@ def test_apply_single_function_pullbacks_triangle3d():
     cov2t = Coefficient(COV2T)
     contra2t = Coefficient(CONTRA2T)
 
+    uml2 = Coefficient(Uml2)
     um = Coefficient(Um)
     vm = Coefficient(Vm)
     vdm = Coefficient(Vdm)
@@ -75,6 +79,7 @@ def test_apply_single_function_pullbacks_triangle3d():
 
     w = Coefficient(W)
 
+    rul2 = ReferenceValue(ul2)
     ru = ReferenceValue(u)
     rv = ReferenceValue(v)
     rvd = ReferenceValue(vd)
@@ -84,6 +89,7 @@ def test_apply_single_function_pullbacks_triangle3d():
     rcov2t = ReferenceValue(cov2t)
     rcontra2t = ReferenceValue(contra2t)
 
+    ruml2 = ReferenceValue(uml2)
     rum = ReferenceValue(um)
     rvm = ReferenceValue(vm)
     rvdm = ReferenceValue(vdm)
@@ -115,6 +121,7 @@ def test_apply_single_function_pullbacks_triangle3d():
 
     mappings = {
         # Simple elements should get a simple representation
+        ul2: rul2 / detJ,
         u: ru,
         v: rv,
         vd: as_vector(M_hdiv[i, j]*rvd[j], i),
@@ -127,6 +134,7 @@ def test_apply_single_function_pullbacks_triangle3d():
         contra2t: as_tensor((1.0 / detJ) * (1.0 / detJ)
                             * J[i, k] * rcontra2t[k, l] * J[j, l], (i, j)),
         # Mixed elements become a bit more complicated
+        uml2: as_vector([ruml2[0] / detJ, ruml2[1] / detJ]),
         um: rum,
         vm: rvm,
         vdm: as_vector([
@@ -204,6 +212,7 @@ def test_apply_single_function_pullbacks_triangle3d():
     }
 
     # Check functions of various elements outside a mixed context
+    check_single_function_pullback(ul2, mappings)
     check_single_function_pullback(u, mappings)
     check_single_function_pullback(v, mappings)
     check_single_function_pullback(vd, mappings)
@@ -214,6 +223,7 @@ def test_apply_single_function_pullbacks_triangle3d():
     check_single_function_pullback(contra2t, mappings)
 
     # Check functions of various elements inside a mixed context
+    check_single_function_pullback(uml2, mappings)
     check_single_function_pullback(um, mappings)
     check_single_function_pullback(vm, mappings)
     check_single_function_pullback(vdm, mappings)
@@ -229,13 +239,15 @@ def test_apply_single_function_pullbacks_triangle():
     cell = triangle
     domain = as_domain(cell)
 
+    Ul2 = FiniteElement("DG L2", cell, 1)
     U = FiniteElement("CG", cell, 1)
     V = VectorElement("CG", cell, 1)
     Vd = FiniteElement("RT", cell, 1)
     Vc = FiniteElement("N1curl", cell, 1)
     T = TensorElement("CG", cell, 1)
     S = TensorElement("CG", cell, 1, symmetry=True)
 
+    Uml2 = Ul2*Ul2
     Um = U*U
     Vm = U*V
     Vdm = V*Vd
@@ -245,13 +257,15 @@ def test_apply_single_function_pullbacks_triangle():
 
     W = S*T*Vc*Vd*V*U
 
+    ul2 = Coefficient(Ul2)
     u = Coefficient(U)
     v = Coefficient(V)
     vd = Coefficient(Vd)
     vc = Coefficient(Vc)
     t = Coefficient(T)
     s = Coefficient(S)
 
+    uml2 = Coefficient(Uml2)
     um = Coefficient(Um)
     vm = Coefficient(Vm)
     vdm = Coefficient(Vdm)
@@ -261,13 +275,15 @@ def test_apply_single_function_pullbacks_triangle():
 
     w = Coefficient(W)
 
+    rul2 = ReferenceValue(ul2)
     ru = ReferenceValue(u)
     rv = ReferenceValue(v)
     rvd = ReferenceValue(vd)
     rvc = ReferenceValue(vc)
     rt = ReferenceValue(t)
     rs = ReferenceValue(s)
 
+    ruml2 = ReferenceValue(uml2)
     rum = ReferenceValue(um)
     rvm = ReferenceValue(vm)
     rvdm = ReferenceValue(vdm)
@@ -294,13 +310,15 @@ def test_apply_single_function_pullbacks_triangle():
 
     mappings = {
         # Simple elements should get a simple representation
+        ul2: rul2 / detJ,
         u: ru,
         v: rv,
         vd: as_vector(M_hdiv[i, j]*rvd[j], i),
         vc: as_vector(Jinv[j, i]*rvc[j], i),
         t: rt,
         s: as_tensor([[rs[0], rs[1]], [rs[1], rs[2]]]),
         # Mixed elements become a bit more complicated
+        uml2: as_vector([ruml2[0] / detJ, ruml2[1] / detJ]),
         um: rum,
         vm: rvm,
         vdm: as_vector([
@@ -358,6 +376,7 @@ def test_apply_single_function_pullbacks_triangle():
     }
 
     # Check functions of various elements outside a mixed context
+    check_single_function_pullback(ul2, mappings)
     check_single_function_pullback(u, mappings)
     check_single_function_pullback(v, mappings)
     check_single_function_pullback(vd, mappings)
@@ -366,6 +385,7 @@ def test_apply_single_function_pullbacks_triangle():
     check_single_function_pullback(s, mappings)
 
     # Check functions of various elements inside a mixed context
+    check_single_function_pullback(uml2, mappings)
     check_single_function_pullback(um, mappings)
     check_single_function_pullback(vm, mappings)
     check_single_function_pullback(vdm, mappings)

test/test_elements.py

@@ -15,7 +15,7 @@
 def test_scalar_galerkin():
     for cell in all_cells:
         for p in range(1, 10):
-            for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG"):
+            for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG", "Discontinuous Lagrange L2", "DG L2"):
                 element = FiniteElement(family, cell, p)
                 assert element.value_shape() == ()
                 assert element == eval(repr(element))
@@ -31,7 +31,7 @@ def test_vector_galerkin():
         # shape = () if dim == 1 else (dim,)
         shape = (dim,)
         for p in range(1, 10):
-            for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG"):
+            for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG", "Discontinuous Lagrange L2", "DG L2"):
                 element = VectorElement(family, cell, p)
                 assert element.value_shape() == shape
                 assert element == eval(repr(element))
@@ -46,7 +46,7 @@ def test_tensor_galerkin():
         # shape = () if dim == 1 else (dim,dim)
         shape = (dim, dim)
         for p in range(1, 10):
-            for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG"):
+            for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG", "Discontinuous Lagrange L2", "DG L2"):
                 element = TensorElement(family, cell, p)
                 assert element.value_shape() == shape
                 assert element == eval(repr(element))
@@ -65,7 +65,7 @@ def test_tensor_symmetry():
                 if isinstance(s, dict) and cell == interval:
                     continue
 
-                for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG"):
+                for family in ("Lagrange", "CG", "Discontinuous Lagrange", "DG", "Discontinuous Lagrange L2", "DG L2"):
                     if isinstance(s, dict):
                         element = TensorElement(
                             family, cell, p, shape=(dim, dim), symmetry=s)
@@ -169,7 +169,8 @@ def test_missing_cell():
         assert element == eval(repr(element))
         element = TensorElement("DG", cell, 1, shape=(2, 2))
         assert element == eval(repr(element))
-
+        element = TensorElement("DG L2", cell, 1, shape=(2, 2))
+        assert element == eval(repr(element))
 
 def test_invalid_degree():
     cell = triangle

ufl/algorithms/apply_function_pullbacks.py

@@ -136,6 +136,9 @@ def apply_single_function_pullbacks(g):
         f = as_tensor((1.0 / detJ) * (1.0 / detJ) * J[i, m] * r[m, n] * J[j, n], (i, j))
         assert f.ufl_shape == g.ufl_shape
         return f
+    elif mapping == "L2 Piola":
+        assert rsh == gsh
+        return r / detJ
 
     # By placing components in a list and using as_vector at the end,
     # we're assuming below that both global function g and its
@@ -225,6 +228,11 @@ def apply_single_function_pullbacks(g):
                                    J[i, m] * rv[m * rdim + n] * J[j, n])
                     g_components[gpos + i * gdim + j] = gv
 
+        elif mp == "L2 Piola":
+            assert gm == rm
+            for i in range(gm):
+                g_components[gpos + i] = r[rpos + i] / detJ
+
         else:
             error("Unknown subelement mapping type %s for element %s." % (mp, str(subelm)))
 

ufl/finiteelement/elementlist.py

@@ -233,6 +233,31 @@ def show_elements():
 register_alias("N2F", lambda family, dim, order,
                degree: ("Brezzi-Douglas-Marini", order))
 
+# discontinuous elements using l2 pullbacks
+register_element2("DPC L2", 0, L2, "L2 Piola", (1, None), cubes)
+register_element2("DQ L2", 0, L2, "L2 Piola", (0, None), cubes)
+register_element("Gauss-Legendre L2", "GL L2", 0, L2, "L2 Piola", (0, None),
+                 ("interval",))
+register_element("Discontinuous Lagrange L2", "DG L2", 0, L2, "L2 Piola", (0, None),
+                 any_cell)  # "DP"
+
+register_alias("DP L2", lambda family, dim, order,
+               degree: ("Discontinuous Lagrange L2", order))
+
+register_alias("P- Lambda L2", lambda family, dim, order,
+               degree: feec_element_l2(family, dim, order, degree))
+register_alias("P Lambda L2", lambda family, dim, order,
+               degree: feec_element_l2(family, dim, order, degree))
+register_alias("Q- Lambda L2", lambda family, dim, order,
+               degree: feec_element_l2(family, dim, order, degree))
+register_alias("S Lambda L2", lambda family, dim, order,
+               degree: feec_element_l2(family, dim, order, degree))
+
+register_alias("P- L2", lambda family, dim, order,
+               degree: feec_element_l2(family, dim, order, degree))
+register_alias("Q- L2", lambda family, dim, order,
+               degree: feec_element_l2(family, dim, order, degree))
+
 
 def feec_element(family, n, r, k):
     """Finite element exterior calculus notation
@@ -280,6 +305,52 @@ def feec_element(family, n, r, k):
     return family, r
 
 
+def feec_element_l2(family, n, r, k):
+    """Finite element exterior calculus notation
+    n = topological dimension of domain
+    r = polynomial order
+    k = form_degree"""
+
+    # Note: We always map to edge elements in 2D, don't know how to
+    # differentiate otherwise?
+
+    # Mapping from (feec name, domain dimension, form degree) to
+    # (family name, polynomial order)
+    _feec_elements = {
+        "P- Lambda L2": (
+            (("P", r), ("DP L2", r - 1)),
+            (("P", r), ("RTE", r), ("DP L2", r - 1)),
+            (("P", r), ("N1E", r), ("N1F", r), ("DP L2", r - 1)),
+        ),
+        "P Lambda L2": (
+            (("P", r), ("DP L2", r)),
+            (("P", r), ("BDME", r), ("DP L2", r)),
+            (("P", r), ("N2E", r), ("N2F", r), ("DP L2", r)),
+        ),
+        "Q- Lambda L2": (
+            (("Q", r), ("DQ L2", r - 1)),
+            (("Q", r), ("RTCE", r), ("DQ L2", r - 1)),
+            (("Q", r), ("NCE", r), ("NCF", r), ("DQ L2", r - 1)),
+        ),
+        "S Lambda L2": (
+            (("S", r), ("DPC L2", r)),
+            (("S", r), ("BDMCE", r), ("DPC L2", r)),
+            (("S", r), ("AAE", r), ("AAF", r), ("DPC L2", r)),
+        ),
+    }
+
+    # New notation, old verbose notation (including "Lambda") might be
+    # removed
+    _feec_elements["P- L2"] = _feec_elements["P- Lambda L2"]
+    _feec_elements["P L2"] = _feec_elements["P Lambda L2"]
+    _feec_elements["Q- L2"] = _feec_elements["Q- Lambda L2"]
+    _feec_elements["S L2"] = _feec_elements["S Lambda L2"]
+
+    family, r = _feec_elements[family][n - 1][k]
+
+    return family, r
+
+
 # General FEEC notation, old verbose (can be removed)
 register_alias("P- Lambda", lambda family, dim, order,
                degree: feec_element(family, dim, order, degree))
@@ -325,6 +396,9 @@ def canonical_element_description(family, cell, order, form_degree):
     if form_degree is not None and family in ("P", "S"):
         family, order = feec_element(family, tdim, order, form_degree)
 
+    if form_degree is not None and family in ("P L2", "S L2"):
+        family, order = feec_element_l2(family, tdim, order, form_degree)
+
     # Check whether this family is an alias for something else
     while family in aliases:
         if tdim is None:
@@ -347,6 +421,10 @@ def canonical_element_description(family, cell, order, form_degree):
             if order >= 1:
                 warning("Discontinuous Lagrange element requested on %s, creating DQ element." % cell.cellname())
             family = "DQ"
+        elif family == "Discontinuous Lagrange L2":
+            if order >= 1:
+                warning("Discontinuous Lagrange L2 element requested on %s, creating DQ L2 element." % cell.cellname())
+            family = "DQ L2"
 
     # Validate cellname if a valid cell is specified
     if not (cellname is None or cellname in cellnames):

ufl/finiteelement/finiteelement.py

@@ -125,6 +125,13 @@ def dq_family(cell):
                                               for c in cell.sub_cells()],
                                             cell=cell)
 
+            elif family == "DQ L2":
+                def dq_family_l2(cell):
+                    return "DG L2" if cell.cellname() in simplices else "DQ L2"
+                return TensorProductElement(*[FiniteElement(dq_family_l2(c), c, degree, variant=variant)
+                                              for c in cell.sub_cells()],
+                                            cell=cell)
+
         return super(FiniteElement, cls).__new__(cls)
 
     def __init__(self,

ufl/finiteelement/tensorproductelement.py

@@ -86,6 +86,8 @@ def __init__(self, *elements, **kwargs):
     def mapping(self):
         if all(e.mapping() == "identity" for e in self._sub_elements):
             return "identity"
+        elif all(e.mapping() == "L2 Piola" for e in self._sub_elements):
+            return "L2 Piola"
         else:
             return "undefined"