Some updates related to dualspace

AUTHORS

@@ -29,3 +29,4 @@ Contributors:
   | Jack S. Hale          <[email protected]>
   | Tuomas Airaksinen     <[email protected]>
   | Reuben W. Hill        <[email protected]>
+  | Nacime Bouziani       <[email protected]>

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, derivative, tetrahedron, triangle, interval, dx)
+                 action, adjoint, derivative, inner, tetrahedron, triangle, interval, dx)
 from ufl.constantvalue import Zero
 from ufl.form import ZeroBaseForm
 
@@ -102,8 +102,6 @@ def test_addition():
     domain_2d = default_domain(triangle)
     f_2d = FiniteElement("CG", triangle, 1)
     V = FunctionSpace(domain_2d, f_2d)
-    f_2d_2 = FiniteElement("CG", triangle, 2)
-    V2 = FunctionSpace(domain_2d, f_2d_2)
     V_dual = V.dual()
 
     u = TrialFunction(V)
@@ -137,11 +135,6 @@ def test_addition():
     res -= ZeroBaseForm((v,))
     assert res == L
 
-    with pytest.raises(ValueError):
-        # Raise error for incompatible arguments
-        v2 = TestFunction(V2)
-        res = L + ZeroBaseForm((v2, u))
-
 
 def test_scalar_mult():
     domain_2d = default_domain(triangle)
@@ -256,7 +249,7 @@ def test_differentiation():
     w = Cofunction(U.dual())
     dwdu = expand_derivatives(derivative(w, u))
     assert isinstance(dwdu, ZeroBaseForm)
-    assert dwdu.arguments() == (Argument(u.ufl_function_space(), 0),)
+    assert dwdu.arguments() == (Argument(w.ufl_function_space().dual(), 0), Argument(u.ufl_function_space(), 1))
     # Check compatibility with int/float
     assert dwdu == 0
 
@@ -285,24 +278,21 @@ def test_differentiation():
     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)
+    Ac = Action(w, u)
+    dAcdu = derivative(Ac, u)
+    assert dAcdu == Action(Adjoint(derivative(w, u)), u) + Action(w, derivative(u, u))
 
-    # -- Adjoint -- #
-    Ad = Adjoint(M)
-    dAddu = expand_derivatives(derivative(Ad, u))
-    # Push differentiation through Adjoint
-    assert dAddu == 0
+    dAcdu = expand_derivatives(dAcdu)
+    # Since dw/du = 0
+    assert dAcdu == 1 * Action(w, v)
 
     # -- Form sum -- #
-    Fs = M + Ac
+    uhat = Argument(U, 1)
+    what = Argument(U, 2)
+    Fs = M + inner(u * uhat, v) * dx
     dFsdu = expand_derivatives(derivative(Fs, u))
     # Distribute differentiation over FormSum components
-    assert dFsdu == 1 * Action(M, v)
+    assert dFsdu == FormSum([inner(what * uhat, v) * dx, 1])
 
 
 def test_zero_base_form_mult():

ufl/action.py

@@ -12,7 +12,8 @@
 from ufl.form import BaseForm, FormSum, Form, ZeroBaseForm
 from ufl.core.ufl_type import ufl_type
 from ufl.algebra import Sum
-from ufl.argument import Argument
+from ufl.constantvalue import Zero
+from ufl.argument import Argument, Coargument
 from ufl.coefficient import BaseCoefficient, Coefficient, Cofunction
 from ufl.differentiation import CoefficientDerivative
 from ufl.matrix import Matrix
@@ -40,19 +41,26 @@ class Action(BaseForm):
         "ufl_operands",
         "_repr",
         "_arguments",
+        "_coefficients",
+        "_domains",
         "_hash")
 
     def __new__(cls, *args, **kw):
         left, right = args
 
         # Check trivial case
         if left == 0 or right == 0:
-            # Check compatibility of function spaces
-            _check_function_spaces(left, right)
             # Still need to work out the ZeroBaseForm arguments.
-            new_arguments = _get_action_form_arguments(left, right)
+            new_arguments, _ = _get_action_form_arguments(left, right)
             return ZeroBaseForm(new_arguments)
 
+        # Coarguments (resp. Argument) from V* to V* (resp. from V to V) are identity matrices,
+        # i.e. we have: V* x V -> R (resp. V x V* -> R).
+        if isinstance(left, (Coargument, Argument)):
+            return right
+        if isinstance(right, (Coargument, Argument)):
+            return left
+
         if isinstance(left, (FormSum, Sum)):
             # Action distributes over sums on the LHS
             return FormSum(*[(Action(component, right), 1)
@@ -70,6 +78,7 @@ def __init__(self, left, right):
         self._left = left
         self._right = right
         self.ufl_operands = (self._left, self._right)
+        self._domains = None
 
         # Check compatibility of function spaces
         _check_function_spaces(left, right)
@@ -97,14 +106,22 @@ def _analyze_form_arguments(self):
         The highest number Argument of the left operand and the lowest number
         Argument of the right operand are consumed by the action.
         """
-        self._arguments = _get_action_form_arguments(self._left, self._right)
+        self._arguments, self._coefficients = _get_action_form_arguments(self._left, self._right)
+
+    def _analyze_domains(self):
+        """Analyze which domains can be found in Action."""
+        from ufl.domain import join_domains
+        # Collect unique domains
+        self._domains = join_domains([e.ufl_domain() for e in self.ufl_operands])
 
     def equals(self, other):
         if type(other) is not Action:
             return False
         if self is other:
             return True
-        return self._left == other._left and self._right == other._right
+        # Make sure we are returning a boolean as left and right equalities can be `ufl.Equation`s
+        # if the underlying objects are `ufl.BaseForm`.
+        return bool(self._left == other._left) and bool(self._right == other._right)
 
     def __str__(self):
         return f"Action({self._left}, {self._right})"
@@ -127,29 +144,51 @@ def _check_function_spaces(left, right):
         # right as a consequence of Leibniz formula.
         right, *_ = right.ufl_operands
 
+    # `left` can also be a Coefficient in V (= V**), e.g. `action(Coefficient(V), Cofunction(V.dual()))`.
+    left_arg = left.arguments()[-1] if not isinstance(left, Coefficient) else left
     if isinstance(right, (Form, Action, Matrix, ZeroBaseForm)):
-        if left.arguments()[-1].ufl_function_space().dual() != right.arguments()[0].ufl_function_space():
+        if left_arg.ufl_function_space().dual() != right.arguments()[0].ufl_function_space():
             raise TypeError("Incompatible function spaces in Action")
     elif isinstance(right, (Coefficient, Cofunction, Argument)):
-        if left.arguments()[-1].ufl_function_space() != right.ufl_function_space():
+        if left_arg.ufl_function_space() != right.ufl_function_space():
             raise TypeError("Incompatible function spaces in Action")
-    else:
+    # `Zero` doesn't contain any information about the function space.
+    # -> Not a problem since Action will get simplified with a `ZeroBaseForm`
+    #    which won't take into account the arguments on the right because of argument contraction.
+    # This occurs for:
+    # `derivative(Action(A, B), u)` with B is an `Expr` such that dB/du == 0
+    # -> `derivative(B, u)` becomes `Zero` when expanding derivatives since B is an Expr.
+    elif not isinstance(right, Zero):
         raise TypeError("Incompatible argument in Action: %s" % type(right))
 
 
 def _get_action_form_arguments(left, right):
     """Perform argument contraction to work out the arguments of Action"""
 
-    if isinstance(right, CoefficientDerivative):
+    coefficients = ()
+    # `left` can also be a Coefficient in V (= V**), e.g. `action(Coefficient(V), Cofunction(V.dual()))`.
+    left_args = left.arguments()[:-1] if not isinstance(left, Coefficient) else ()
+    if isinstance(right, BaseForm):
+        arguments = left_args + right.arguments()[1:]
+        coefficients += right.coefficients()
+    elif isinstance(right, CoefficientDerivative):
         # Action differentiation pushes differentiation through
         # right as a consequence of Leibniz formula.
-        right, *_ = right.ufl_operands
-
-    if isinstance(right, BaseForm):
-        return left.arguments()[:-1] + right.arguments()[1:]
-    elif isinstance(right, BaseCoefficient):
-        return left.arguments()[:-1]
+        from ufl.algorithms.analysis import extract_arguments_and_coefficients
+        right_args, right_coeffs = extract_arguments_and_coefficients(right)
+        arguments = left_args + tuple(right_args)
+        coefficients += tuple(right_coeffs)
+    elif isinstance(right, (BaseCoefficient, Zero)):
+        arguments = left_args
+        # When right is ufl.Zero, Action gets simplified so updating
+        # coefficients here doesn't matter
+        coefficients += (right,)
     elif isinstance(right, Argument):
-        return left.arguments()[:-1] + (right,)
+        arguments = left_args + (right,)
     else:
         raise TypeError
+
+    if isinstance(left, BaseForm):
+        coefficients += left.coefficients()
+
+    return arguments, coefficients

ufl/adjoint.py

@@ -10,6 +10,7 @@
 # Modified by Nacime Bouziani, 2021-2022.
 
 from ufl.form import BaseForm, FormSum, ZeroBaseForm
+from ufl.argument import Coargument
 from ufl.core.ufl_type import ufl_type
 # --- The Adjoint class represents the adjoint of a numerical object that
 #     needs to be computed at assembly time ---
@@ -28,6 +29,8 @@ class Adjoint(BaseForm):
         "_form",
         "_repr",
         "_arguments",
+        "_coefficients",
+        "_domains",
         "ufl_operands",
         "_hash")
 
@@ -44,6 +47,14 @@ def __new__(cls, *args, **kw):
             # Adjoint distributes over sums
             return FormSum(*[(Adjoint(component), 1)
                              for component in form.components()])
+        elif isinstance(form, Coargument):
+            # The adjoint of a coargument `c: V* -> V*` is the identity matrix mapping from V to V (i.e. V x V* -> R).
+            # Equivalently, the adjoint of `c` is its first argument, which is a ufl.Argument defined on the
+            # primal space of `c`.
+            primal_arg, _ = form.arguments()
+            # Returning the primal argument avoids explicit argument reconstruction, making it
+            # a robust strategy for handling subclasses of `ufl.Coargument`.
+            return primal_arg
 
         return super(Adjoint, cls).__new__(cls)
 
@@ -55,6 +66,7 @@ def __init__(self, form):
 
         self._form = form
         self.ufl_operands = (self._form,)
+        self._domains = None
         self._hash = None
         self._repr = "Adjoint(%s)" % repr(self._form)
 
@@ -68,13 +80,22 @@ def form(self):
     def _analyze_form_arguments(self):
         """The arguments of adjoint are the reverse of the form arguments."""
         self._arguments = self._form.arguments()[::-1]
+        self._coefficients = self._form.coefficients()
+
+    def _analyze_domains(self):
+        """Analyze which domains can be found in Adjoint."""
+        from ufl.domain import join_domains
+        # Collect unique domains
+        self._domains = join_domains([e.ufl_domain() for e in self.ufl_operands])
 
     def equals(self, other):
         if type(other) is not Adjoint:
             return False
         if self is other:
             return True
-        return self._form == other._form
+        # Make sure we are returning a boolean as the equality can result in a `ufl.Equation`
+        # if the underlying objects are `ufl.BaseForm`.
+        return bool(self._form == other._form)
 
     def __str__(self):
         return f"Adjoint({self._form})"

ufl/algorithms/__init__.py

@@ -21,6 +21,7 @@
     "estimate_total_polynomial_degree",
     "sort_elements",
     "compute_form_data",
+    "preprocess_form",
     "apply_transformer",
     "ReuseTransformer",
     "load_ufl_file",
@@ -79,7 +80,7 @@
 
 # Preprocessing a form to extract various meta data
 # from ufl.algorithms.formdata import FormData
-from ufl.algorithms.compute_form_data import compute_form_data
+from ufl.algorithms.compute_form_data import compute_form_data, preprocess_form
 
 # Utilities for checking properties of forms
 from ufl.algorithms.signature import compute_form_signature

ufl/algorithms/ad.py

@@ -11,7 +11,6 @@
 
 import warnings
 
-from ufl.adjoint import Adjoint
 from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering
 from ufl.algorithms.apply_derivatives import apply_derivatives
 
@@ -29,13 +28,6 @@ def expand_derivatives(form, **kwargs):
     if kwargs:
         warnings("Deprecation: expand_derivatives no longer takes any keyword arguments")
 
-    if isinstance(form, Adjoint):
-        dform = expand_derivatives(form._form)
-        if dform == 0:
-            return dform
-        # Adjoint is taken on a 3-form which can't happen
-        raise NotImplementedError('Adjoint derivative is not supported.')
-
     # Lower abstractions for tensor-algebra types into index notation
     form = apply_algebra_lowering(form)
 

ufl/algorithms/analysis.py

@@ -50,12 +50,14 @@ def extract_type(a, ufl_types):
     if not isinstance(ufl_types, (list, tuple)):
         ufl_types = (ufl_types,)
 
-    # BaseForms that aren't forms only have arguments
+    # BaseForms that aren't forms only contain arguments & coefficients
     if isinstance(a, BaseForm) and not isinstance(a, Form):
+        objects = set()
         if any(issubclass(t, BaseArgument) for t in ufl_types):
-            return set(a.arguments())
-        else:
-            return set()
+            objects.update(a.arguments())
+        if any(issubclass(t, BaseCoefficient) for t in ufl_types):
+            objects.update(a.coefficients())
+        return objects
 
     if all(issubclass(t, Terminal) for t in ufl_types):
         # Optimization

ufl/algorithms/apply_derivatives.py

@@ -24,7 +24,7 @@
 from ufl.core.terminal import Terminal
 from ufl.corealg.map_dag import map_expr_dag
 from ufl.corealg.multifunction import MultiFunction
-from ufl.differentiation import CoordinateDerivative
+from ufl.differentiation import CoordinateDerivative, BaseFormCoordinateDerivative
 from ufl.domain import extract_unique_domain
 from ufl.operators import (bessel_I, bessel_J, bessel_K, bessel_Y, cell_avg,
                            conditional, cos, cosh, exp, facet_avg, ln, sign,
@@ -1032,7 +1032,7 @@ def cofunction(self, o):
         dc = self.coefficient(o)
         if dc == 0:
             # Convert ufl.Zero into ZeroBaseForm
-            return ZeroBaseForm(self._v)
+            return ZeroBaseForm(o.arguments() + self._v)
         return dc
 
     def coargument(self, o):
@@ -1102,6 +1102,14 @@ def coordinate_derivative(self, o, f, dummy_w, dummy_v, dummy_cd):
                                                  rcache=self.rcaches[key]),
                                     o_[1], o_[2], o_[3])
 
+    def base_form_coordinate_derivative(self, o, f, dummy_w, dummy_v, dummy_cd):
+        o_ = o.ufl_operands
+        key = (BaseFormCoordinateDerivative, o_[0])
+        return BaseFormCoordinateDerivative(map_expr_dag(self, o_[0],
+                                                         vcache=self.vcaches[key],
+                                                         rcache=self.rcaches[key]),
+                                            o_[1], o_[2], o_[3])
+
     def indexed(self, o, Ap, ii):  # TODO: (Partially) duplicated in generic rules
         # Reuse if untouched
         if Ap is o.ufl_operands[0]: