From b3068a72b6f1356570cf16e0e8bcb0a74368bed2 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 18 Jun 2026 21:20:14 +0100 Subject: [PATCH 1/2] FormProduct --- test/test_form.py | 107 +++++++++++++++++++++++++++ ufl/__init__.py | 3 +- ufl/algorithms/map_integrands.py | 9 ++- ufl/algorithms/traversal.py | 5 +- ufl/classes.py | 3 +- ufl/form.py | 120 ++++++++++++++++++++++++++++++- ufl/formoperators.py | 23 +++++- 7 files changed, 263 insertions(+), 7 deletions(-) diff --git a/test/test_form.py b/test/test_form.py index 0d680f82e..a777d808e 100755 --- a/test/test_form.py +++ b/test/test_form.py @@ -2,9 +2,11 @@ from utils import LagrangeElement from ufl import ( + Argument, Coefficient, Cofunction, Form, + FormProduct, FormSum, FunctionSpace, Mesh, @@ -19,6 +21,7 @@ grad, inner, nabla_grad, + replace, triangle, ) from ufl.form import BaseForm @@ -203,3 +206,107 @@ def test_formsum(mass): assert f.weights()[0] == -1 assert isinstance(df, FormSum) assert df.weights()[0] == -9 + + +def test_form_product_constructor_and_arguments(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + v = TestFunction(V) + f = Coefficient(V) + g = Coefficient(V) + h = Coefficient(V) + + Lf = f * v * dx + Lg = g * v * dx + Lh = h * v * dx + + product = FormProduct(Lf, Lg) + assert isinstance(product, BaseForm) + assert product.factors() == (Lf, Lg) + assert product.ufl_operands == (Lf, Lg) + assert product.factor_arguments() == (Lf.arguments(), Lg.arguments()) + + arguments = product.arguments() + assert tuple(argument.number() for argument in arguments) == (0, 1) + assert tuple(argument.part() for argument in arguments) == (None, None) + assert tuple(argument.ufl_function_space() for argument in arguments) == (V, V) + assert Lg.arguments()[0].number() == 0 + + assert product.coefficients() == (f, g) + assert product.ufl_domains() == (domain,) + + nested = FormProduct(Lf, FormProduct(Lg, Lh)) + assert nested.factors() == (Lf, Lg, Lh) + assert tuple(argument.number() for argument in nested.arguments()) == (0, 1, 2) + + +def test_form_product_rejects_invalid_inputs(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + v = TestFunction(V) + f = Coefficient(V) + L = f * v * dx + + with pytest.raises(ValueError): + FormProduct(L) + with pytest.raises(TypeError): + FormProduct(L, 1) + + +def test_form_product_is_explicit_not_mul_overload(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + v = TestFunction(V) + f = Coefficient(V) + g = Coefficient(V) + Lf = f * v * dx + Lg = g * v * dx + + with pytest.raises(TypeError): + Lf * Lg + + +def test_form_product_replace(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + v = TestFunction(V) + f = Coefficient(V) + g = Coefficient(V) + + Lf = f * v * dx + Lg = g * v * dx + product = FormProduct(Lf, Lg) + replaced = replace(product, {f: g}) + + assert isinstance(replaced, FormProduct) + assert bool(replaced.factors()[0] == Lg) + assert bool(replaced.factors()[1] == Lg) + assert tuple(argument.number() for argument in replaced.arguments()) == (0, 1) + + +def test_form_product_derivative_product_rule(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + v = TestFunction(V) + f = Coefficient(V) + direction = Argument(V, 1) + + L = f * v * dx + product = FormProduct(L, L) + dproduct = derivative(product, f, direction) + + assert isinstance(dproduct, FormSum) + assert len(dproduct.components()) == 2 + assert all(isinstance(component, FormProduct) for component in dproduct.components()) + + dL = derivative(L, f, direction) + expected_first = FormProduct(dL, L) + expected_second = FormProduct(L, dL) + assert bool(dproduct.components()[0] == expected_first) + assert bool(dproduct.components()[1] == expected_second) + + +def test_form_product_exported_from_classes(): + from ufl.classes import FormProduct as ClassesFormProduct + + assert ClassesFormProduct is FormProduct diff --git a/ufl/__init__.py b/ufl/__init__.py index 415914d77..41ecb70eb 100644 --- a/ufl/__init__.py +++ b/ufl/__init__.py @@ -271,7 +271,7 @@ from ufl.core.multiindex import Index, indices from ufl.domain import AbstractDomain, Mesh, MeshSequence, MeshView from ufl.finiteelement import AbstractFiniteElement -from ufl.form import BaseForm, Form, FormSum, ZeroBaseForm +from ufl.form import BaseForm, Form, FormProduct, FormSum, ZeroBaseForm from ufl.formoperators import ( action, adjoint, @@ -491,6 +491,7 @@ "FacetArea", "FacetNormal", "Form", + "FormProduct", "FormSum", "FunctionSpace", "H1Curl", diff --git a/ufl/algorithms/map_integrands.py b/ufl/algorithms/map_integrands.py index ddf70979b..a810e77ab 100644 --- a/ufl/algorithms/map_integrands.py +++ b/ufl/algorithms/map_integrands.py @@ -15,7 +15,7 @@ from ufl.constantvalue import Zero from ufl.core.expr import Expr from ufl.corealg.map_dag import map_expr_dag -from ufl.form import BaseForm, Form, FormSum, ZeroBaseForm +from ufl.form import BaseForm, Form, FormProduct, FormSum, ZeroBaseForm from ufl.integral import Integral @@ -69,6 +69,13 @@ def map_integrands(function, form, only_integral_type=None): right = map_integrands(function, form._right, only_integral_type) # Zeros are caught inside `Action.__new__` return Action(left, right) + elif isinstance(form, FormProduct): + factors = tuple( + map_integrands(function, factor, only_integral_type) for factor in form.factors() + ) + if any(factor == 0 for factor in factors): + return ZeroBaseForm(form.arguments()) + return FormProduct(*factors) elif isinstance(form, ZeroBaseForm): arguments = tuple( map_integrands(function, arg, only_integral_type) for arg in form._arguments diff --git a/ufl/algorithms/traversal.py b/ufl/algorithms/traversal.py index 8fcc84239..c6d7661a0 100644 --- a/ufl/algorithms/traversal.py +++ b/ufl/algorithms/traversal.py @@ -11,7 +11,7 @@ from ufl.action import Action from ufl.adjoint import Adjoint from ufl.core.expr import Expr -from ufl.form import BaseForm, Form, FormSum +from ufl.form import BaseForm, Form, FormProduct, FormSum from ufl.integral import Integral @@ -23,6 +23,7 @@ def iter_expressions(a): - a is an Integral: the integrand expression of a - a is a Form: all integrand expressions of all integrals - a is a FormSum: the components of a + - a is a FormProduct: the factors of a - a is an Action: the left and right component of a - a is an Adjoint: the underlying form of a """ @@ -30,7 +31,7 @@ def iter_expressions(a): return (itg.integrand() for itg in a.integrals()) elif isinstance(a, Integral): return (a.integrand(),) - elif isinstance(a, FormSum | Adjoint | Action): + elif isinstance(a, FormSum | FormProduct | Adjoint | Action): return tuple(e for op in a.ufl_operands for e in iter_expressions(op)) elif isinstance(a, Expr | BaseForm): return (a,) diff --git a/ufl/classes.py b/ufl/classes.py index bf8463d55..6db866917 100644 --- a/ufl/classes.py +++ b/ufl/classes.py @@ -88,7 +88,7 @@ from ufl.equation import Equation from ufl.exprcontainers import ExprList, ExprMapping from ufl.finiteelement import AbstractFiniteElement -from ufl.form import BaseForm, Form, FormSum, ZeroBaseForm +from ufl.form import BaseForm, Form, FormProduct, FormSum, ZeroBaseForm from ufl.functionspace import ( AbstractFunctionSpace, DualSpace, @@ -323,6 +323,7 @@ "Form", "Form", "FormArgument", + "FormProduct", "FormSum", "FunctionSpace", "GeometricCellQuantity", diff --git a/ufl/form.py b/ufl/form.py index 9388d8f25..7a3cc8090 100644 --- a/ufl/form.py +++ b/ufl/form.py @@ -33,7 +33,7 @@ from ufl.classes import AbstractDomain # Export list for ufl.classes -__all_classes__ = ["Form", "BaseForm", "ZeroBaseForm"] +__all_classes__ = ["Form", "BaseForm", "FormProduct", "ZeroBaseForm"] # --- The Form class, representing a complete variational form or functional --- @@ -697,6 +697,124 @@ def as_form(form): return form +@ufl_type() +class FormProduct(BaseForm): + """Form product. + + A structural product of variational forms and form-like objects. + Factors are stored unchanged; only the product's aggregate argument + list is renumbered cumulatively across factor argument lists. + """ + + __slots__ = ( + "_arguments", + "_coefficients", + "_domains", + "_factors", + "_geometric_quantities", + "_hash", + "ufl_operands", + ) + _ufl_required_methods_ = "_analyze_form_arguments" # type: ignore + + def __init__(self, *factors): + """Initialise.""" + BaseForm.__init__(self) + + if len(factors) < 2: + raise ValueError("FormProduct requires at least two factors.") + + full_factors = [] + for factor in factors: + if isinstance(factor, FormProduct): + full_factors.extend(factor.factors()) + elif isinstance(factor, BaseForm): + full_factors.append(factor) + else: + raise TypeError(f"Expected a UFL BaseForm instance, got {type(factor)}.") + + self._factors = tuple(full_factors) + self.ufl_operands = self._factors + self._arguments = None + self._coefficients = None + self._geometric_quantities = None + self._domains = None + self._hash = None + + def factors(self): + """Return the unchanged product factors.""" + return self._factors + + def factor_arguments(self): + """Return original factor-local arguments for each factor.""" + return tuple(factor.arguments() for factor in self._factors) + + def _analyze_form_arguments(self): + """Analyze aggregate arguments and coefficients.""" + arguments = [] + coefficients = [] + cumulative_shift = 0 + + for factor in self._factors: + factor_arguments = factor.arguments() + arguments.extend( + type(argument)( + argument.ufl_function_space(), + argument.number() + cumulative_shift, + argument.part(), + ) + for argument in factor_arguments + ) + coefficients.extend(factor.coefficients()) + cumulative_shift += len(factor_arguments) + + self._arguments = tuple(arguments) + self._coefficients = tuple(sorted(set(coefficients), key=lambda x: x.count())) + self._geometric_quantities = () + + def _analyze_domains(self): + """Analyze which domains can be found in FormProduct.""" + from ufl.domain import join_domains, sort_domains + + self._domains = sort_domains( + join_domains(chain.from_iterable(factor.ufl_domains() for factor in self._factors)) + ) + + def ufl_domains(self): + """Return all domains found in the base form.""" + if self._domains is None: + self._analyze_domains() + return self._domains + + def equals(self, other): + """Evaluate ``bool(lhs_form_product == rhs_form_product)``.""" + if type(other) is not FormProduct: + return False + if self is other: + return True + return len(self.factors()) == len(other.factors()) and all( + bool(a == b) for a, b in zip(self.factors(), other.factors()) + ) + + def __hash__(self): + """Hash.""" + if self._hash is None: + self._hash = hash(("FormProduct", tuple(hash(factor) for factor in self.factors()))) + return self._hash + + def empty(self): + """Returns whether any product factor is empty.""" + return any(factor.empty() for factor in self.factors()) + + def __str__(self): + """Compute shorter string representation of form product.""" + return "FormProduct({})".format(", ".join(str(factor) for factor in self._factors)) + + def __repr__(self): + """Representation.""" + return "FormProduct({})".format(", ".join(repr(factor) for factor in self._factors)) + + @ufl_type() class FormSum(BaseForm): """Form sum. diff --git a/ufl/formoperators.py b/ufl/formoperators.py index 59d6ded47..2bca76764 100644 --- a/ufl/formoperators.py +++ b/ufl/formoperators.py @@ -42,7 +42,7 @@ ) from ufl.exprcontainers import ExprList, ExprMapping from ufl.finiteelement import AbstractFiniteElement -from ufl.form import BaseForm, Form, FormSum, ZeroBaseForm, as_form +from ufl.form import BaseForm, Form, FormProduct, FormSum, ZeroBaseForm, as_form from ufl.functionspace import FunctionSpace from ufl.geometry import SpatialCoordinate from ufl.indexed import Indexed @@ -397,6 +397,27 @@ def derivative(form, coefficient, argument=None, coefficient_derivatives=None): for component, w in zip(form.components(), form.weights()) ] ) + elif isinstance(form, FormProduct): + # Apply the product rule while reusing one derivative argument across all factors. + _, arguments = _handle_derivative_arguments(form, coefficient, argument) + derivative_argument = tuple(arguments) + product_terms = [] + factors = form.factors() + for i, factor in enumerate(factors): + differentiated_factor = derivative( + factor, coefficient, derivative_argument, coefficient_derivatives + ) + if differentiated_factor == 0 or isinstance(differentiated_factor, ZeroBaseForm): + continue + if isinstance(differentiated_factor, BaseForm) and differentiated_factor.empty(): + continue + term_factors = list(factors) + term_factors[i] = differentiated_factor + product_terms.append((FormProduct(*term_factors), 1)) + + if not product_terms: + return ZeroBaseForm(form.arguments() + derivative_argument) + return FormSum(*product_terms) elif isinstance(form, Adjoint): # Is `derivative(Adjoint(A), ...)` with A a 2-form even legal ? # -> If yes, what's the right thing to do here ? From addea580fea8fe31e22271a9f35fe2e57d964aa1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 19 Jun 2026 05:24:43 +0100 Subject: [PATCH 2/2] Adjoint(FormProduct) --- test/test_form.py | 48 ++++++++++++++++++++++++++++++++++++++++++++++- ufl/adjoint.py | 10 +++++++++- ufl/form.py | 9 +++++++++ 3 files changed, 65 insertions(+), 2 deletions(-) diff --git a/test/test_form.py b/test/test_form.py index a777d808e..c3cf31d7c 100755 --- a/test/test_form.py +++ b/test/test_form.py @@ -2,6 +2,7 @@ from utils import LagrangeElement from ufl import ( + Adjoint, Argument, Coefficient, Cofunction, @@ -9,6 +10,7 @@ FormProduct, FormSum, FunctionSpace, + Matrix, Mesh, SpatialCoordinate, TestFunction, @@ -240,6 +242,16 @@ def test_form_product_constructor_and_arguments(domain): assert tuple(argument.number() for argument in nested.arguments()) == (0, 1, 2) +def test_form_product_of_one_factor_simplifies(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + v = TestFunction(V) + f = Coefficient(V) + L = f * v * dx + + assert FormProduct(L) is L + + def test_form_product_rejects_invalid_inputs(domain): element = LagrangeElement(triangle, 1) V = FunctionSpace(domain, element) @@ -248,7 +260,9 @@ def test_form_product_rejects_invalid_inputs(domain): L = f * v * dx with pytest.raises(ValueError): - FormProduct(L) + FormProduct() + with pytest.raises(TypeError): + FormProduct(1) with pytest.raises(TypeError): FormProduct(L, 1) @@ -266,6 +280,38 @@ def test_form_product_is_explicit_not_mul_overload(domain): Lf * Lg +def test_adjoint_form_product_reverses_adjoint_factors(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + A = Matrix(V, V) + B = Matrix(V, V) + C = Matrix(V, V) + + product = FormProduct(A, B, C) + adjoint_product = Adjoint(product) + + assert isinstance(adjoint_product, FormProduct) + assert tuple(factor.form() for factor in adjoint_product.factors()) == (C, B, A) + assert adjoint_product.factors() == (Adjoint(C), Adjoint(B), Adjoint(A)) + + +def test_adjoint_form_product_leaves_rank_zero_and_one_factors_unadjointed(domain): + element = LagrangeElement(triangle, 1) + V = FunctionSpace(domain, element) + v = TestFunction(V) + f = Coefficient(V) + functional = f * dx + linear = f * v * dx + A = Matrix(V, V) + + product = FormProduct(functional, linear, A) + adjoint_product = Adjoint(product) + + assert isinstance(adjoint_product, FormProduct) + assert adjoint_product.factors() == (Adjoint(A), linear, functional) + assert Adjoint(FormProduct(functional, linear)).factors() == (linear, functional) + + def test_form_product_replace(domain): element = LagrangeElement(triangle, 1) V = FunctionSpace(domain, element) diff --git a/ufl/adjoint.py b/ufl/adjoint.py index d4c170cb6..63ba5b7a4 100644 --- a/ufl/adjoint.py +++ b/ufl/adjoint.py @@ -12,7 +12,7 @@ from ufl.argument import Coargument from ufl.core.ufl_type import ufl_type -from ufl.form import BaseForm, FormSum, ZeroBaseForm +from ufl.form import BaseForm, FormProduct, FormSum, ZeroBaseForm # --- The Adjoint class represents the adjoint of a numerical object that # needs to be computed at assembly time --- @@ -50,6 +50,14 @@ def __new__(cls, *args, **kw): elif isinstance(form, FormSum): # Adjoint distributes over sums return FormSum(*((Adjoint(c), w) for c, w in zip(form.components(), form.weights()))) + elif isinstance(form, FormProduct): + # Reverse product order and take the adjoint of rank-2 factors. + return FormProduct( + *( + factor if len(factor.arguments()) < 2 else Adjoint(factor) + for factor in reversed(form.factors()) + ) + ) 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). diff --git a/ufl/form.py b/ufl/form.py index 7a3cc8090..7bf6000ab 100644 --- a/ufl/form.py +++ b/ufl/form.py @@ -717,6 +717,15 @@ class FormProduct(BaseForm): ) _ufl_required_methods_ = "_analyze_form_arguments" # type: ignore + def __new__(cls, *factors): + """Create a new FormProduct.""" + if len(factors) == 1: + (factor,) = factors + if isinstance(factor, BaseForm): + return factor + raise TypeError(f"Expected a UFL BaseForm instance, got {type(factor)}.") + return super().__new__(cls) + def __init__(self, *factors): """Initialise.""" BaseForm.__init__(self)