Dualspace update

test/test_duals.py

@@ -4,7 +4,7 @@
 from ufl import FiniteElement, FunctionSpace, MixedFunctionSpace, \
     Coefficient, Matrix, Cofunction, FormSum, Argument, Coargument,\
     TestFunction, TrialFunction, Adjoint, Action, \
-    action, adjoint, tetrahedron, triangle, interval, dx
+    action, adjoint, derivative, tetrahedron, triangle, interval, dx
 
 __authors__ = "India Marsden"
 __date__ = "2020-12-28 -- 2020-12-28"
@@ -13,6 +13,7 @@
 
 from ufl.domain import default_domain
 from ufl.duals import is_primal, is_dual
+from ufl.algorithms.ad import expand_derivatives
 
 
 def test_mixed_functionspace(self):
@@ -58,7 +59,6 @@ def test_dual_coefficients():
     v = Coefficient(V, count=1)
     u = Coefficient(V_dual, count=1)
     w = Cofunction(V_dual)
-    x = Cofunction(V)
 
     assert is_primal(v)
     assert not is_dual(v)
@@ -69,8 +69,11 @@ def test_dual_coefficients():
     assert is_dual(w)
     assert not is_primal(w)
 
-    assert is_primal(x)
-    assert not is_dual(x)
+    try:
+        x = Cofunction(V)
+        assert False
+    except ValueError:
+        pass
 
 
 def test_dual_arguments():
@@ -82,7 +85,6 @@ def test_dual_arguments():
     v = Argument(V, 1)
     u = Argument(V_dual, 2)
     w = Coargument(V_dual, 3)
-    x = Coargument(V, 4)
 
     assert is_primal(v)
     assert not is_dual(v)
@@ -93,8 +95,11 @@ def test_dual_arguments():
     assert is_dual(w)
     assert not is_primal(w)
 
-    assert is_primal(x)
-    assert not is_dual(x)
+    try:
+        x = Coargument(V, 4)
+        assert False
+    except ValueError:
+        pass
 
 
 def test_addition():
@@ -153,6 +158,9 @@ def test_adjoint():
     assert isinstance(res, FormSum)
     assert isinstance(res.components()[0], Adjoint)
 
+    # Adjoint(Adjoint(.)) = Id
+    assert adjoint(adj) == a
+
 
 def test_action():
     domain_2d = default_domain(triangle)
@@ -167,7 +175,7 @@ def test_action():
     u = Coefficient(U)
     u_a = Argument(U, 0)
     v = Coefficient(V)
-    u_star = Cofunction(U.dual())
+    ustar = Cofunction(U.dual())
     u_form = u_a * dx
 
     res = action(a, u)
@@ -191,4 +199,57 @@ def test_action():
         res = action(a, v)
 
     with pytest.raises(TypeError):
-        res = action(a, u_star)
+        res = action(a, ustar)
+
+
+def test_differentiation():
+    domain_2d = default_domain(triangle)
+    f_2d = FiniteElement("CG", triangle, 1)
+    V = FunctionSpace(domain_2d, f_2d)
+    domain_1d = default_domain(interval)
+    f_1d = FiniteElement("CG", interval, 1)
+    U = FunctionSpace(domain_1d, f_1d)
+
+    u = Coefficient(U)
+    v = Argument(U, 0)
+    vstar = Argument(U.dual(), 0)
+
+    # -- Cofunction -- #
+    w = Cofunction(U.dual())
+    dwdu = expand_derivatives(derivative(w, u))
+    assert dwdu == 0
+
+    dwdw = expand_derivatives(derivative(w, w, vstar))
+    assert dwdw == vstar
+
+    dudw = expand_derivatives(derivative(u, w))
+    assert dudw == 0
+
+    # -- Coargument -- #
+    dvstardu = expand_derivatives(derivative(vstar, u))
+    assert dvstardu == 0
+
+    # -- Matrix -- #
+    M = Matrix(V, U)
+    dMdu = expand_derivatives(derivative(M, u))
+    assert dMdu == 0
+
+    # -- Action -- #
+    Ac = Action(M, u)
+    dAcdu = expand_derivatives(derivative(Ac, u))
+
+    # Action(dM/du, u) + Action(M, du/du) = Action(M, uhat) since dM/du = 0.
+    # Multiply by 1 to get a FormSum (type compatibility).
+    assert dAcdu == 1 * Action(M, v)
+
+    # -- Adjoint -- #
+    Ad = Adjoint(M)
+    dAddu = expand_derivatives(derivative(Ad, u))
+    # Push differentiation through Adjoint
+    assert dAddu == 0
+
+    # -- Form sum -- #
+    Fs = M + Ac
+    dFsdu = expand_derivatives(derivative(Fs, u))
+    # Distribute differentiation over FormSum components
+    assert dFsdu == 1 * Action(M, v)

ufl/__init__.py

@@ -344,7 +344,7 @@
 from ufl.measure import Measure, register_integral_type, integral_types, custom_integral_types
 
 # Form class
-from ufl.form import Form, FormSum, replace_integral_domains
+from ufl.form import Form, BaseForm, FormSum, replace_integral_domains
 
 # Integral classes
 from ufl.integral import Integral

ufl/action.py

@@ -8,7 +8,9 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 from ufl.form import BaseForm, FormSum, Form
-from ufl.coefficient import Coefficient
+from ufl.argument import Argument
+from ufl.coefficient import Coefficient, Cofunction
+from ufl.differentiation import CoefficientDerivative
 
 # --- The Action class represents the action of a numerical object that needs
 #     to be computed at assembly time ---
@@ -29,6 +31,7 @@ class Action(BaseForm):
     __slots__ = (
         "_left",
         "_right",
+        "ufl_operands",
         "_repr",
         "_arguments",
         "_hash")
@@ -39,6 +42,10 @@ def __getnewargs__(self):
     def __new__(cls, *args, **kw):
         left, right = args
 
+        # Check trivial case
+        if left == 0 or right == 0:
+            return 0
+
         if isinstance(left, FormSum):
             # Action distributes over sums on the LHS
             return FormSum(*[(Action(component, right), 1)
@@ -55,19 +62,25 @@ def __init__(self, left, right):
 
         self._left = left
         self._right = right
+        self.ufl_operands = (self._left, self._right)
+
+        if isinstance(right, CoefficientDerivative):
+            # Action differentiation pushes differentiation through
+            # right as a consequence of Leibniz formula.
+            right, *_ = right.ufl_operands
 
-        if isinstance(right, Form):
+        if isinstance(right, (Form, Action)):
             if (left.arguments()[-1].ufl_function_space().dual()
                 != right.arguments()[0].ufl_function_space()):
 
                 raise TypeError("Incompatible function spaces in Action")
-        elif isinstance(right, Coefficient):
+        elif isinstance(right, (Coefficient, Cofunction, Argument)):
             if (left.arguments()[-1].ufl_function_space()
                 != right.ufl_function_space()):
 
                 raise TypeError("Incompatible function spaces in Action")
         else:
-            raise TypeError("Incompatible argument in Action")
+            raise TypeError("Incompatible argument in Action: %s" % type(right))
 
         self._repr = "Action(%s, %s)" % (repr(self._left), repr(self._right))
         self._hash = None
@@ -101,6 +114,13 @@ def _analyze_form_arguments(self):
         else:
             raise TypeError
 
+    def __eq__(self, other):
+        if not isinstance(other, Action):
+            return False
+        if self is other:
+            return True
+        return (self._left == other._left and self._right == other._right)
+
     def __str__(self):
         return "Action(%s, %s)" % (repr(self._left), repr(self._right))
 

ufl/adjoint.py

@@ -23,14 +23,21 @@ class Adjoint(BaseForm):
         "_form",
         "_repr",
         "_arguments",
+        "ufl_operands",
         "_hash")
 
     def __getnewargs__(self):
         return (self._form)
 
     def __new__(cls, *args, **kw):
         form = args[0]
-        if isinstance(form, FormSum):
+        # Check trivial case
+        if form == 0:
+            return 0
+
+        if isinstance(form, Adjoint):
+            return form._form
+        elif isinstance(form, FormSum):
             # Adjoint distributes over sums
             return FormSum(*[(Adjoint(component), 1)
                              for component in form.components()])
@@ -44,6 +51,7 @@ def __init__(self, form):
             raise ValueError("Can only take Adjoint of a 2-form.")
 
         self._form = form
+        self.ufl_operands = (self._form,)
         self._hash = None
         self._repr = "Adjoint(%s)" % repr(self._form)
 
@@ -58,6 +66,13 @@ def _analyze_form_arguments(self):
         """The arguments of adjoint are the reverse of the form arguments."""
         self._arguments = self._form.arguments()[::-1]
 
+    def __eq__(self, other):
+        if not isinstance(other, Adjoint):
+            return False
+        if self is other:
+            return True
+        return (self._form == other._form)
+
     def __str__(self):
         return "Adjoint(%s)" % self._form
 

ufl/algorithms/ad.py

@@ -10,6 +10,9 @@
 # Modified by Anders Logg, 2009.
 
 from ufl.log import warning
+from ufl.form import FormSum
+from ufl.action import Action
+from ufl.adjoint import Adjoint
 from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering
 from ufl.algorithms.apply_derivatives import apply_derivatives
 
@@ -27,6 +30,17 @@ def expand_derivatives(form, **kwargs):
     if kwargs:
         warning("Deprecation: expand_derivatives no longer takes any keyword arguments")
 
+    if isinstance(form, FormSum):
+        return FormSum(*[(expand_derivatives(component), 1) for component in form.components()])
+    if isinstance(form, Action):
+        return Action(expand_derivatives(form._left), expand_derivatives(form._right))
+    if isinstance(form, Adjoint):
+        dform = expand_derivatives(form._form)
+        if not dform:
+            # If dform == 0
+            return dform
+        raise NotImplementedError('Adjoint derivative is not supported.')
+
     # Lower abstractions for tensor-algebra types into index notation
     form = apply_algebra_lowering(form)