schnellerhase · Sep 8, 2026 · Oct 7, 2026 · Oct 7, 2026 · Oct 8, 2026 · Oct 12, 2026
diff --git a/ufl/algorithms/analysis.py b/ufl/algorithms/analysis.py
@@ -22,6 +22,7 @@
 from ufl.domain import Mesh
 from ufl.form import BaseForm, Form
 from ufl.geometry import GeometricQuantity
+from ufl.integral import Integral
 from ufl.utils.sorting import sorted_by_count, topological_sorting
 
 # TODO: Some of these can possibly be optimised by implementing
@@ -71,7 +72,11 @@ def extract_type(a, ufl_types, base_form_op_as_expr=None):
         remove_base_form_ops = False
 
     if base_form_op_as_expr is None:
-        base_form_op_as_expr = isinstance(a, Form | Expr) and not isinstance(a, BaseFormOperator)
+        # An integral is part of a form, so a base form operator in one has its
+        # Coargument bound just as it does in the form.
+        base_form_op_as_expr = isinstance(a, Form | Integral | Expr) and not isinstance(
+            a, BaseFormOperator
+        )
 
     # BaseForms that aren't forms or base form operators
     # only contain arguments & coefficients

diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py
@@ -955,8 +955,9 @@ def _(self, o: Expr) -> Expr:
 
         Represent grad(grad(f)) as Grad(Grad(f)).
         """
-        # Check that o is a "differential terminal"
-        if not isinstance(o.ufl_operands[0], Grad | Terminal):
+        # Check that o is a "differential terminal". A base form operator is
+        # one too: the handler above leaves `grad(N)` unevaluated.
+        if not isinstance(o.ufl_operands[0], Grad | Terminal | BaseFormOperator):
             raise ValueError("Expecting only grads applied to a terminal.")
         return Grad(o)
 
@@ -1274,11 +1275,13 @@ def __init__(
         compress: bool | None = True,
         visited_cache: dict[tuple, Expr | BaseForm] | None = None,
         result_cache: dict[Expr | BaseForm, Expr | BaseForm] | None = None,
+        differentiate_through_operators: bool = False,
     ) -> None:
         """Initialise."""
         super().__init__(
             (), compress=compress, visited_cache=visited_cache, result_cache=result_cache
         )
+        self._differentiate_through_operators = differentiate_through_operators
         # Type checking
         if not isinstance(coefficients, ExprList):
             raise ValueError("Expecting a ExprList of coefficients.")
@@ -1593,6 +1596,20 @@ def compute_gprimeterm(ngrads, vval, vcomp, wshape, wcomp):
                     else:
                         gprimesum += prod
 
+        if isinstance(o, BaseFormOperator) and isinstance(gprimesum, Zero):
+            # `grad(N)` is differentiated as if `N` were a Coefficient, so the
+            # partial derivative above is zero. The chain rule through `N` must
+            # not be dropped: either push the derivative into the operator, or
+            # record it so that `apply_derivatives` adds the
+            # Action(dF/dN, dN/dw) term.
+            if self._differentiate_through_operators:
+                do = self(o)
+                # A zero carries no domain, so it must not be wrapped in `Grad`.
+                if isinstance(do, Zero):
+                    return Zero(g.ufl_shape)
+                return apply_grads(do)
+            self.pending_operations += (o,)
+
         return gprimesum
 
     @process.register(CoordinateDerivative)
@@ -1614,8 +1631,22 @@ def _(self, o: BaseFormOperator, *dfs) -> Expr:
         """
         d_coeff = self._process_coefficient(o)
         # It also handles the non-scalar case
-        if d_coeff == 0:
-            self.pending_operations += (o,)
+        if d_coeff != 0:
+            return d_coeff
+
+        target_space = getattr(o, "target_space", None)
+        if self._differentiate_through_operators and target_space is not None:
+            # A base form operator is linear in its operand, so the Gateaux
+            # derivative commutes with it: D_w[v](N(g(w))) = N(D_w[v](g(w))).
+            # Rebuilding on the target space rather than the existing Coargument
+            # lets the operator number the Coargument away from the new
+            # argument, which higher derivatives need.
+            (df,) = dfs
+            if isinstance(df, Zero):
+                return Zero(o.ufl_shape)
+            return o._ufl_expr_reconstruct_(df, v=target_space())
+
+        self.pending_operations += (o,)
         return d_coeff
 
     # -- Handlers for BaseForm objects -- #
@@ -1780,9 +1811,11 @@ def __init__(
         compress: bool | None = True,
         visited_cache: dict[tuple, Expr | BaseForm] | None = None,
         result_cache: dict[Expr | BaseForm, Expr | BaseForm] | None = None,
+        differentiate_through_operators: bool = False,
     ) -> None:
         """Initialise."""
         super().__init__(compress=compress, visited_cache=visited_cache, result_cache=result_cache)
+        self._differentiate_through_operators = differentiate_through_operators
         # Record the operations delayed to the derivative expansion phase:
         # Example: dN(u)/du where `N` is a BaseFormOperator and `u` a Coefficient
         self.pending_operations = ()
@@ -1860,7 +1893,12 @@ def _(self, o: CoefficientDerivative, f: Expr | BaseForm) -> Expr | BaseForm:
         # operations
         dag_traverser = self._dag_traverser_cache.setdefault(
             key,
-            GateauxDerivativeRuleset(w, v, cd),  # type: ignore
+            GateauxDerivativeRuleset(  # type: ignore
+                w,
+                v,
+                cd,  # type: ignore
+                differentiate_through_operators=self._differentiate_through_operators,
+            ),
         )
         # If f has been seen by the traverser, it immediately returns
         # the cached value.
@@ -1996,18 +2034,26 @@ def __iadd__(self, other):
         return self
 
 
-def apply_derivatives(expression):
+def apply_derivatives(expression, differentiate_through_operators: bool = False):
     """Apply derivatives to an expression.
 
     Args:
         expression: A Form, an Expr or a BaseFormOperator to be differentiated
+        differentiate_through_operators: If True, differentiate a base form
+            operator by differentiating its operand, which is exact because the
+            operator is linear in it, instead of recording the derivative as an
+            action to be assembled separately. This composes to higher
+            derivatives and keeps the result an ordinary form, which a form
+            compiler that evaluates the operator cell-locally needs.
 
     Returns:
         A differentiated expression
     """
     # Notation: Let `var` be the thing we are differentating with respect to.
 
-    dag_traverser = DerivativeRuleDispatcher()
+    dag_traverser = DerivativeRuleDispatcher(
+        differentiate_through_operators=differentiate_through_operators
+    )
 
     # If we hit a base form operator (bfo), then if `var` is:
     #    - a BaseFormOperator → Return `d(expression)/dw` where `w` is

diff --git a/ufl/algorithms/balancing.py b/ufl/algorithms/balancing.py
@@ -9,6 +9,7 @@
 from ufl.classes import (
     Grad,
     Indexed,
+    Interpolate,
     NegativeRestricted,
     PositiveRestricted,
     ReferenceGrad,
@@ -32,11 +33,21 @@
 }
 
 
+def _is_modifiable(expr):
+    """Check whether modifiers may be applied to ``expr``.
+
+    An interpolation is modifiable like a terminal: it is a finite element field
+    of its target space, and the expression it interpolates is evaluated
+    separately.
+    """
+    return expr._ufl_is_terminal_ or isinstance(expr, Interpolate)
+
+
 def balance_modified_terminal(expr):
     """Balance modified terminal."""
     # NB! Assuming e.g. grad(cell_avg(expr)) does not occur,
     # i.e. it is simplified to 0 immediately.
-    if expr._ufl_is_terminal_:
+    if _is_modifiable(expr):
         return expr
 
     assert expr._ufl_is_terminal_modifier_
@@ -45,12 +56,12 @@ def balance_modified_terminal(expr):
 
     # Build list of modifier layers
     layers = [expr]
-    while not expr._ufl_is_terminal_:
+    while not _is_modifiable(expr):
         assert expr._ufl_is_terminal_modifier_
         expr = expr.ufl_operands[0]
         layers.append(expr)
     assert layers[-1] is expr
-    assert expr._ufl_is_terminal_
+    assert _is_modifiable(expr)
 
     # Apply modifiers in order
     layers = sorted(layers[:-1], key=lambda e: modifier_precedence[e._ufl_handler_name_])

diff --git a/ufl/core/interpolate.py b/ufl/core/interpolate.py
@@ -55,7 +55,9 @@ def __init__(self, expr, v):
         if isinstance(v, AbstractFunctionSpace):
             if is_dual(v):
                 raise ValueError("Expecting a primal function space.")
-            n = 1 if expr_arg_numbers == {0} else 0
+            # Take the lowest free number, so that an expression linear in
+            # several arguments stays representable.
+            n = min(set(range(len(expr_arg_numbers) + 1)) - expr_arg_numbers)
             v = Argument(v.dual(), n)
             dual_arg_numbers = {n}
         elif isinstance(v, BaseForm):
@@ -66,7 +68,8 @@ def __init__(self, expr, v):
         # Check valid argument numbering
         if expr_arg_numbers & dual_arg_numbers:
             raise ValueError("Same argument numbers in first and second operands to interpolate.")
-        if expr_arg_numbers | dual_arg_numbers not in [set(), {0}, {0, 1}]:
+        argument_numbers = expr_arg_numbers | dual_arg_numbers
+        if argument_numbers and argument_numbers != set(range(len(argument_numbers))):
             raise ValueError("Non-contiguous argument numbers in interpolate.")
 
         # Reversed order convention
@@ -155,6 +158,15 @@ def signature(slot):
             self._signature = hashlib.sha512(str(signatures).encode("utf-8")).hexdigest()
         return self._signature
 
+    def target_space(self) -> AbstractFunctionSpace:
+        """Return the function space that the expression is interpolated into.
+
+        Unlike `ufl_function_space` this does not depend on the argument
+        numbering, so it stays correct when the interpolated expression holds an
+        argument of the enclosing form.
+        """
+        return self._function_space
+
     def ufl_element(self) -> AbstractFiniteElement:
         """Return the target finite element."""
         return self._function_space.ufl_element()

diff --git a/ufl/differentiation.py b/ufl/differentiation.py
@@ -19,6 +19,7 @@
 from ufl.exprcontainers import ExprList, ExprMapping
 from ufl.form import BaseForm
 from ufl.precedence import parstr
+from ufl.utils.counted import Counted
 from ufl.variable import Variable
 
 # --- Basic differentiation objects ---
@@ -157,6 +158,11 @@ def __init__(self, base_form, coefficients, arguments, coefficient_derivatives):
         )
         self._argument_slots = base_form._argument_slots
         self._domains = None
+        # `BaseFormOperator` is a `Counted`, but the `BaseFormDerivative.__init__`
+        # above does not reach `Counted.__init__`, so initialise the count here.
+        # Without it `Form.terminal_numbering` (and hence `Form.signature`) fails
+        # on any form holding a derivative of a base form operator.
+        Counted.__init__(self, counted_class=BaseFormOperator)
 
     # Enforce Operator reconstruction as Operator is a parent class of
     # both: BaseFormDerivative and BaseFormOperator.