Merge master (Constant changes) into mixed-dimensional + fix conflits

doc/sphinx/source/manual/form_language.rst

@@ -424,9 +424,8 @@ There is a shorthand for this, whose use is similar to ``Arguments``, called
 
   u, p = Coefficients(TH)
 
-Spatially constant (or discontinuous piecewise constant) functions can
-conveniently be represented by ``Constant``, ``VectorConstant``, and
-``TensorConstant``::
+Spatially constant values can conveniently be represented by
+``Constant``, ``VectorConstant``, and ``TensorConstant``::
 
   c0 = Constant(cell)
   v0 = VectorConstant(cell)

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_derivative.py

@@ -603,18 +603,15 @@ def Nw(x, derivatives):
         assert derivatives == (0,)
         return dw
 
-    w, b, K = form_data_f.original_form.coefficients()
     mapping = {K: Kv, b: bv, w: Nw}
     fv2 = f_expression((0,), mapping)
     self.assertAlmostEqual(fv, fv2)
 
-    w, b, K = form_data_F.original_form.coefficients()
     v, = form_data_F.original_form.arguments()
     mapping = {K: Kv, b: bv, v: Nv, w: Nw}
     Fv2 = F_expression((0,), mapping)
     self.assertAlmostEqual(Fv, Fv2)
 
-    w, b, K = form_data_J.original_form.coefficients()
     v, u = form_data_J.original_form.arguments()
     mapping = {K: Kv, b: bv, v: Nv, u: Nu, w: Nw}
     Jv2 = J_expression((0,), mapping)

test/test_diff.py

@@ -173,7 +173,10 @@ def df(v):
 
 
 def testCoefficient():
-    v = Constant(triangle)
+    coord_elem = VectorElement("P", triangle, 1, dim=3)
+    mesh = Mesh(coord_elem)
+    V = FunctionSpace(mesh, FiniteElement("P", triangle, 1))
+    v = Coefficient(V)
     assert round(expand_derivatives(diff(v, v))-1.0, 7) == 0
 
 

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

@@ -301,8 +301,8 @@
     Arguments, TestFunctions, TrialFunctions
 
 # Coefficients
-from ufl.coefficient import Coefficient, Coefficients, \
-    Constant, VectorConstant, TensorConstant
+from ufl.coefficient import Coefficient, Coefficients
+from ufl.constant import Constant, VectorConstant, TensorConstant
 
 # Split function
 from ufl.split_functions import split

ufl/algorithms/analysis.py

@@ -29,6 +29,7 @@
 from ufl.core.terminal import Terminal, FormArgument
 from ufl.argument import Argument
 from ufl.coefficient import Coefficient
+from ufl.constant import Constant
 from ufl.algorithms.traversal import iter_expressions
 from ufl.corealg.traversal import unique_pre_traversal, traverse_unique_terminals
 
@@ -110,6 +111,11 @@ def extract_coefficients(a):
     return sorted_by_count(extract_type(a, Coefficient))
 
 
+def extract_constants(a):
+    """Build a sorted list of all constants in a"""
+    return sorted_by_count(extract_type(a, Constant))
+
+
 def extract_arguments_and_coefficients(a):
     """Build two sorted lists of all arguments and coefficients
     in a, which can be a Form, Integral or Expr."""

ufl/algorithms/apply_derivatives.py

@@ -116,6 +116,9 @@ def independent_operator(self, o):
     # Literals are by definition independent of any differentiation variable
     constant_value = independent_terminal
 
+    # Constants are independent of any differentiation
+    constant = independent_terminal
+
     # Rules for form arguments must be specified in specialized rule set
     form_argument = override
 

ufl/algorithms/apply_function_pullbacks.py

@@ -33,6 +33,7 @@
 
 from ufl.tensors import as_tensor, as_vector
 from ufl.utils.sequences import product
+import numpy
 
 
 def sub_elements_with_mappings(element):
@@ -135,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
@@ -146,14 +150,12 @@ def apply_single_function_pullbacks(g):
     #     (ONLY IF REFERENCE VALUE SHAPE PRESERVES TENSOR RANK)
     #   - All cases with scalar subelements and without symmetries are
     #     covered by the shortcut above
-    # - VectorElements of vector-valued basic elements (FIXME)
-    # - TensorElements with symmetries (FIXME)
-    assert len(gsh) == 1
-    assert len(rsh) == 1
 
     g_components = [None] * gsize
     gpos = 0
     rpos = 0
+
+    r = as_vector([r[idx] for idx in numpy.ndindex(r.ufl_shape)])
     for subelm in sub_elements_with_mappings(element):
         gm = product(subelm.value_shape())
         rm = product(subelm.reference_value_shape())
@@ -226,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)))
 
@@ -234,8 +241,7 @@ def apply_single_function_pullbacks(g):
 
     # Wrap up components in a vector, must return same shape as input
     # function g
-    assert len(gsh) == 1
-    f = as_vector(g_components)
+    f = as_tensor(numpy.asarray(g_components).reshape(gsh))
     assert f.ufl_shape == g.ufl_shape
     return f
 

ufl/algorithms/estimate_degrees.py

@@ -50,6 +50,9 @@ def constant_value(self, v):
         "Constant values are constant."
         return 0
 
+    def constant(self, v):
+        return 0
+
     def geometric_quantity(self, v):
         "Some geometric quantities are cellwise constant. Others are nonpolynomial and thus hard to estimate."
         if is_cellwise_constant(v):

ufl/algorithms/formfiles.py

@@ -30,6 +30,7 @@
 from ufl.form import Form
 from ufl.finiteelement import FiniteElementBase
 from ufl.core.expr import Expr
+from ufl.constant import Constant
 from ufl.argument import Argument
 from ufl.coefficient import Coefficient
 
@@ -148,7 +149,7 @@ def interpret_ufl_namespace(namespace):
             # FIXME: Remove after FFC is updated to use reserved_objects:
             ufd.object_names[name] = value
             ufd.object_by_name[name] = value
-        elif isinstance(value, (FiniteElementBase, Coefficient, Argument, Form, Expr)):
+        elif isinstance(value, (FiniteElementBase, Coefficient, Constant, Argument, Form, Expr)):
             # Store instance <-> name mappings for important objects
             # without a reserved name
             ufd.object_names[id(value)] = name

ufl/algorithms/signature.py

@@ -22,7 +22,7 @@
 from ufl.classes import (Label,
                          Index, MultiIndex,
                          Coefficient, Argument,
-                         GeometricQuantity, ConstantValue,
+                         GeometricQuantity, ConstantValue, Constant,
                          ExprList, ExprMapping)
 from ufl.log import error
 from ufl.corealg.traversal import traverse_unique_terminals, pre_traversal
@@ -71,6 +71,9 @@ def compute_terminal_hashdata(expressions, renumbering):
             elif isinstance(expr, Coefficient):
                 data = expr._ufl_signature_data_(renumbering)
 
+            elif isinstance(expr, Constant):
+                data = expr._ufl_signature_data_(renumbering)
+
             elif isinstance(expr, Argument):
                 data = expr._ufl_signature_data_(renumbering)