From a48aa01168c3f2daa4ab9640d590533b9460d352 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 13:52:50 +0100 Subject: [PATCH 01/35] initial class --- ufl/domain.py | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/ufl/domain.py b/ufl/domain.py index 9ba4372e2..c92a7b5cb 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -10,7 +10,7 @@ import numbers from collections.abc import Iterable, Sequence -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Union if TYPE_CHECKING: from ufl.core.expr import Expr @@ -21,6 +21,7 @@ from ufl.core.ufl_type import UFLObject from ufl.corealg.traversal import traverse_unique_terminals from ufl.sobolevspace import H1 +from ufl.corealg.dag_traverser import DAGTraverser # Export list for ufl.classes __all_classes__ = ["AbstractDomain", "Mesh", "MeshView"] @@ -444,3 +445,21 @@ def find_geometric_dimension(expr): raise ValueError("Cannot determine geometric dimension from expression.") (gdim,) = gdims return gdim + + +class UniqueDomainExtractor(DAGTraverser): + # Terminals: Arg or Coeff or constant. constant has no domain. + # Operators: 1) Children match domain, return that. 2) Differing domains, raise error. 3) Different domains, return A or B + + # ExternalOperator or Interpolate: Domain is where we are mapping into. These are both BaseFormOperators. + def __init__( + self, + compress: Union[bool, None] = True, + visited_cache: Union[dict[tuple, Expr], None] = None, + result_cache: Union[dict[Expr, Expr], None] = None, + ) -> None: + """Initialise.""" + self._compress = compress + self._visited_cache = {} if visited_cache is None else visited_cache + self._result_cache = {} if result_cache is None else result_cache + super().__init__(compress=compress, visited_cache=visited_cache, result_cache=result_cache) From 13711c89f5c1a44bdae82d074604de9e164497af Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 13:56:50 +0100 Subject: [PATCH 02/35] add test --- test/test_unique_domain_extractor.py | 49 ++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 test/test_unique_domain_extractor.py diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py new file mode 100644 index 000000000..44ab275aa --- /dev/null +++ b/test/test_unique_domain_extractor.py @@ -0,0 +1,49 @@ +from utils import FiniteElement, LagrangeElement, MixedElement + +import pytest +from ufl import ( + Coefficient, + FunctionSpace, + Mesh, + MeshSequence, + SpatialCoordinate, + TrialFunction, + split, + triangle, + cos +) + +from ufl.domain import extract_unique_domain_dag +from ufl.pullback import contravariant_piola, identity_pullback +from ufl.sobolevspace import L2, HDiv + +def test_extract_unique_domain(): + cell = triangle + elem0 = LagrangeElement(cell, 1) + elem1 = FiniteElement("Brezzi-Douglas-Marini", cell, 2, (2,), contravariant_piola, HDiv) + elem2 = FiniteElement("Discontinuous Lagrange", cell, 1, (), identity_pullback, L2) + elem = MixedElement([elem0, elem1, elem2]) + mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=100) + mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=101) + mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=102) + domain = MeshSequence([mesh1, mesh2, mesh3]) + V = FunctionSpace(domain, elem) + + u = TrialFunction(V) + u1, u2, u3 = split(u) + for i, u_i in enumerate((u1, u2, u3)): + assert extract_unique_domain_dag(u_i) == domain[i] + + f = Coefficient(V) + f1, f2, f3 = split(f) + for i, f_i in enumerate((f1, f2, f3)): + assert extract_unique_domain_dag(f_i) == domain[i] + + x1, y1 = SpatialCoordinate(mesh1) + expr = u1 + x1 * cos(x1) + assert extract_unique_domain_dag(expr) == mesh1 + + x2, y2 = SpatialCoordinate(mesh2) + with pytest.raises(ValueError) as e_info: + _ = extract_unique_domain_dag(u1 + u2) + _ = extract_unique_domain_dag(u1 + u2 + x2 * cos(x2 * u1)) \ No newline at end of file From 14e250e6345a40e6bab80211a4a3a21e3a3fc252 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 14:57:09 +0100 Subject: [PATCH 03/35] add function --- ufl/domain.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/ufl/domain.py b/ufl/domain.py index c92a7b5cb..665f41288 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -463,3 +463,39 @@ def __init__( self._visited_cache = {} if visited_cache is None else visited_cache self._result_cache = {} if result_cache is None else result_cache super().__init__(compress=compress, visited_cache=visited_cache, result_cache=result_cache) + +def extract_unique_domain_dag(expr: Union[Expr, Form]) -> AbstractDomain: + """Extract the single unique domain from an expression. + + This works for expressions containing Indexed Arguments and Coefficients from + split functions on mixed function spaces. + + Args: + expr: Expr or Form to extract domain from + + Returns: + AbstractDomain: The unique domain extracted from the expression. + """ + from ufl.form import Form + from ufl.core.expr import Expr + + if isinstance(expr, Form): + # For forms, we extract domains from integrals + domains = set() + for integral in expr.integrals(): + domain = extract_unique_domain_dag(integral.integrand()) + if domain is not None: + domains.add(domain) + + if len(domains) == 0: + return None + elif len(domains) == 1: + return domains[0] + else: + raise ValueError(f"Form has multiple domains: {domains}") + elif isinstance(expr, Expr): + # For expressions, use the DAG traverser + extractor = UniqueDomainExtractor() + return extractor(expr) + else: + raise TypeError(f"Unsupported type for domain extraction: {type(expr)}. Expected Expr or Form.") \ No newline at end of file From 2cf60e403dfc63b1db0ddc4f33efc4e9da6dff27 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 15:01:14 +0100 Subject: [PATCH 04/35] working for simple expressions --- ufl/checks.py | 2 +- ufl/domain.py | 95 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 96 insertions(+), 1 deletion(-) diff --git a/ufl/checks.py b/ufl/checks.py index 132bb7559..9dab24615 100644 --- a/ufl/checks.py +++ b/ufl/checks.py @@ -11,7 +11,6 @@ from ufl.core.expr import Expr from ufl.core.terminal import FormArgument from ufl.corealg.traversal import traverse_unique_terminals -from ufl.geometry import GeometricQuantity from ufl.sobolevspace import H1 @@ -40,6 +39,7 @@ def is_cellwise_constant(expr): def is_scalar_constant_expression(expr): """Check if an expression is a globally constant scalar expression.""" + from ufl.geometry import GeometricQuantity if is_python_scalar(expr): return True if expr.ufl_shape: diff --git a/ufl/domain.py b/ufl/domain.py index 665f41288..00a1cf979 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -11,6 +11,7 @@ import numbers from collections.abc import Iterable, Sequence from typing import TYPE_CHECKING, Union +from functools import singledispatchmethod if TYPE_CHECKING: from ufl.core.expr import Expr @@ -22,6 +23,9 @@ from ufl.corealg.traversal import traverse_unique_terminals from ufl.sobolevspace import H1 from ufl.corealg.dag_traverser import DAGTraverser +from ufl.core.operator import Operator +from ufl.core.terminal import Terminal +from ufl.indexed import Indexed # Export list for ufl.classes __all_classes__ = ["AbstractDomain", "Mesh", "MeshView"] @@ -464,6 +468,97 @@ def __init__( self._result_cache = {} if result_cache is None else result_cache super().__init__(compress=compress, visited_cache=visited_cache, result_cache=result_cache) + @singledispatchmethod + def process(self, o: Expr) -> Expr: + """Process ``o``. + + Args: + o: `Expr` to be processed. + + Returns: + Processed object. + + """ + return super().process(o) + + @process.register(Indexed) + @DAGTraverser.postorder + def _(self, o: Expr, *operand_results) -> AbstractDomain: + """Process Indexed object by extracting the domain corresponding to the index.""" + from ufl.functionspace import FunctionSpace + # Get the indexed expression and the multiindex + expression, multiindex = o.ufl_operands + + # Get the domain from the first operand (the expression being indexed) + expression_domain = operand_results[0] + + # If the expression domain is a MeshSequence, extract the specific mesh + if isinstance(expression_domain, MeshSequence): + index = multiindex[0]._value + # For mixed function spaces with MeshSequence, we need to determine + # which component mesh this index corresponds to + # This requires understanding the mixed element structure + element = expression.ufl_element() + if hasattr(element, 'sub_elements'): + # Calculate which sub-element (and thus which mesh) this index belongs to + offset = 0 + for i, sub_element in enumerate(element.sub_elements): + # Get the value size for this sub-element on its corresponding mesh + mesh_for_element = expression_domain.meshes[i] + sub_element_fs = FunctionSpace(mesh_for_element, sub_element) + sub_element_size = sub_element_fs.value_size + + if index < offset + sub_element_size: + # This index belongs to mesh i + return mesh_for_element + offset += sub_element_size + raise ValueError(f"Index {index} out of range for mixed function space") + else: + # Simple case: direct indexing into MeshSequence + if 0 <= index < len(expression_domain.meshes): + return expression_domain.meshes[index] + + # If it's not a MeshSequence, just return the expression domain + return expression_domain + + @process.register(Terminal) + def _(self, o: Expr) -> AbstractDomain: + from ufl.coefficient import Coefficient + from ufl.argument import Argument + if isinstance(o, (Coefficient, Argument)): + fs = o.ufl_function_space() + return fs.ufl_domain() + try: + return o.ufl_domain() + except AttributeError: + # If the terminal does not have a domain, return None + return None + + @process.register(Operator) + @DAGTraverser.postorder + def _(self, o: Expr, *operand_results) -> AbstractDomain: + """Process Operator.""" + # Filter out None results (from operands that don't have domains) + domains = [d for d in operand_results if d is not None] + + if not domains: + # No operands have domains + return None + elif len(domains) == 1: + # Only one operand has a domain + return domains[0] + else: + # Multiple operands have domains - they should all be the same + first_domain = domains[0] + if all(d == first_domain for d in domains): + return first_domain + else: + # If domains are not none and differ, raise error + raise ValueError( + f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" + ) + + def extract_unique_domain_dag(expr: Union[Expr, Form]) -> AbstractDomain: """Extract the single unique domain from an expression. From 8937403e32cf0d4c0d1ec20650cacbff39df55fb Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 15:02:23 +0100 Subject: [PATCH 05/35] fix --- ufl/domain.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ufl/domain.py b/ufl/domain.py index 00a1cf979..b4225e1f6 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -585,7 +585,7 @@ def extract_unique_domain_dag(expr: Union[Expr, Form]) -> AbstractDomain: if len(domains) == 0: return None elif len(domains) == 1: - return domains[0] + return domains.pop() else: raise ValueError(f"Form has multiple domains: {domains}") elif isinstance(expr, Expr): From d8d2f6f2e89d63097a1fdb2ff2f61cbff8b00635 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 15:02:31 +0100 Subject: [PATCH 06/35] add form test --- test/test_unique_domain_extractor.py | 31 ++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index 44ab275aa..215a57e87 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -7,10 +7,12 @@ Mesh, MeshSequence, SpatialCoordinate, + Measure, TrialFunction, split, triangle, - cos + cos, + inner ) from ufl.domain import extract_unique_domain_dag @@ -46,4 +48,29 @@ def test_extract_unique_domain(): x2, y2 = SpatialCoordinate(mesh2) with pytest.raises(ValueError) as e_info: _ = extract_unique_domain_dag(u1 + u2) - _ = extract_unique_domain_dag(u1 + u2 + x2 * cos(x2 * u1)) \ No newline at end of file + _ = extract_unique_domain_dag(u1 + u2 + x2 * cos(x2 * u1)) + + +def test_extract_unique_domain_form(): + cell = triangle + elem0 = LagrangeElement(cell, 1) + elem1 = FiniteElement("Brezzi-Douglas-Marini", cell, 2, (2,), contravariant_piola, HDiv) + elem2 = FiniteElement("Discontinuous Lagrange", cell, 1, (), identity_pullback, L2) + elem = MixedElement([elem0, elem1, elem2]) + mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=100) + mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=101) + mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=102) + domain = MeshSequence([mesh1, mesh2, mesh3]) + V = FunctionSpace(domain, elem) + + u = TrialFunction(V) + u1, u2, u3 = split(u) + + f = Coefficient(V) + f1, f2, f3 = split(f) + + dx1 = Measure("dx", mesh1) + + form1 = inner(u1, f1) * dx1 + + assert extract_unique_domain_dag(form1) == mesh1 \ No newline at end of file From 5b20d6df184b7c549c553997a0de77583acabc97 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 15:04:41 +0100 Subject: [PATCH 07/35] add more tests --- test/test_unique_domain_extractor.py | 192 ++++++++++++++++++++++++++- 1 file changed, 190 insertions(+), 2 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index 215a57e87..384a47604 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -12,7 +12,8 @@ split, triangle, cos, - inner + inner, + div ) from ufl.domain import extract_unique_domain_dag @@ -73,4 +74,191 @@ def test_extract_unique_domain_form(): form1 = inner(u1, f1) * dx1 - assert extract_unique_domain_dag(form1) == mesh1 \ No newline at end of file + assert extract_unique_domain_dag(form1) == mesh1 + + +def test_extract_unique_domain_single_mesh(): + """Test domain extraction for standard function spaces on a single mesh.""" + cell = triangle + mesh = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=200) + + # Test scalar elements + P1 = LagrangeElement(cell, 1) + V_scalar = FunctionSpace(mesh, P1) + u_scalar = TrialFunction(V_scalar) + f_scalar = Coefficient(V_scalar) + + # Test that single mesh functions return the correct domain + assert extract_unique_domain_dag(u_scalar) == mesh + assert extract_unique_domain_dag(f_scalar) == mesh + + # Test vector elements + P1_vec = LagrangeElement(cell, 1, (2,)) + V_vector = FunctionSpace(mesh, P1_vec) + u_vector = TrialFunction(V_vector) + f_vector = Coefficient(V_vector) + + assert extract_unique_domain_dag(u_vector) == mesh + assert extract_unique_domain_dag(f_vector) == mesh + + # Test indexing into vector elements + assert extract_unique_domain_dag(u_vector[0]) == mesh + assert extract_unique_domain_dag(u_vector[1]) == mesh + assert extract_unique_domain_dag(f_vector[0]) == mesh + assert extract_unique_domain_dag(f_vector[1]) == mesh + + # Test tensor elements + P1_tensor = LagrangeElement(cell, 1, (2, 2)) + V_tensor = FunctionSpace(mesh, P1_tensor) + u_tensor = TrialFunction(V_tensor) + f_tensor = Coefficient(V_tensor) + + assert extract_unique_domain_dag(u_tensor) == mesh + assert extract_unique_domain_dag(f_tensor) == mesh + assert extract_unique_domain_dag(u_tensor[0, 0]) == mesh + assert extract_unique_domain_dag(u_tensor[1, 1]) == mesh + assert extract_unique_domain_dag(f_tensor[0, 1]) == mesh + + # Test expressions combining same-domain functions + x, y = SpatialCoordinate(mesh) + expr1 = u_scalar + f_scalar + expr2 = u_vector[0] + x + expr3 = inner(u_vector, f_vector) + + assert extract_unique_domain_dag(expr1) == mesh + assert extract_unique_domain_dag(expr2) == mesh + assert extract_unique_domain_dag(expr3) == mesh + + # Test forms + dx = Measure("dx", mesh) + form = inner(u_scalar, f_scalar) * dx + assert extract_unique_domain_dag(form) == mesh + + +def test_extract_unique_domain_mixed_scalar_vector_tensor(): + """Test domain extraction for mixed function spaces with scalar, vector, and tensor elements.""" + cell = triangle + mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) + mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) + mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) + domain = MeshSequence([mesh1, mesh2, mesh3]) + + # Create mixed element with different types: + # - scalar on mesh1 + # - vector on mesh2 + # - tensor on mesh3 + scalar_elem = LagrangeElement(cell, 1) # 1 component + vector_elem = LagrangeElement(cell, 1, (2,)) # 2 components + tensor_elem = LagrangeElement(cell, 1, (2, 2)) # 4 components + mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) + + V = FunctionSpace(domain, mixed_elem) + u = TrialFunction(V) + f = Coefficient(V) + + # For mixed element [scalar, vector(2), tensor(2,2)], split gives: + # - u_components[0]: scalar (mesh1) - index 0 + # - u_components[1]: vector (mesh2) - indices 1-2 + # - u_components[2]: tensor (mesh3) - indices 3-6 + u_scalar, u_vector, u_tensor = split(u) + f_scalar, f_vector, f_tensor = split(f) + + # Test that each component maps to correct mesh + for i, u_i in enumerate((u_scalar, u_vector, u_tensor)): + assert extract_unique_domain_dag(u_i) == domain[i] + for i, f_i in enumerate((f_scalar, f_vector, f_tensor)): + assert extract_unique_domain_dag(f_i) == domain[i] + + # Test indexing into vector and tensor components + for i in range(2): + assert extract_unique_domain_dag(u_vector[i]) == mesh2 + assert extract_unique_domain_dag(f_vector[i]) == mesh2 + + for i in range(2): + for j in range(2): + assert extract_unique_domain_dag(u_tensor[i, j]) == mesh3 + assert extract_unique_domain_dag(f_tensor[i, j]) == mesh3 + + # Test expressions on same mesh (should work) + x1, y1 = SpatialCoordinate(mesh1) + x2, y2 = SpatialCoordinate(mesh2) + x3, y3 = SpatialCoordinate(mesh3) + + # Scalar expressions + expr_scalar = u_scalar * y1 + f_scalar + x1 + assert extract_unique_domain_dag(expr_scalar) == mesh1 + + # Vector expressions + expr_vector = inner(u_vector * y2, f_vector) + x2 + assert extract_unique_domain_dag(expr_vector) == mesh2 + + # Vector component expressions + expr_vec_comp = u_vector[0] + f_vector[1] * y2 + x2 + assert extract_unique_domain_dag(expr_vec_comp) == mesh2 + + # Tensor expressions + expr_tensor = y3 * u_tensor[0, 0] + f_tensor[1, 1] + x3 + assert extract_unique_domain_dag(expr_tensor) == mesh3 + + # Test expressions mixing different mesh components (should fail) + with pytest.raises(ValueError, match="Cannot extract unique domain"): + extract_unique_domain_dag(u_scalar + u_vector[0]) + + with pytest.raises(ValueError, match="Cannot extract unique domain"): + extract_unique_domain_dag(u_vector[0] + u_tensor[0, 0]) + + with pytest.raises(ValueError, match="Cannot extract unique domain"): + extract_unique_domain_dag(f_scalar + f_tensor[1, 1]) + + with pytest.raises(ValueError, match="Cannot extract unique domain"): + extract_unique_domain_dag(u_scalar + x2) + + with pytest.raises(ValueError, match="Cannot extract unique domain"): + extract_unique_domain_dag(u_vector[0] + x3) + + # Test forms on individual meshes + dx1 = Measure("dx", mesh1) + dx2 = Measure("dx", mesh2) + dx3 = Measure("dx", mesh3) + + form_scalar = u_scalar * f_scalar * dx1 + form_vector = inner(u_vector, f_vector) * dx2 + form_tensor = u_tensor[0, 0] * f_tensor[1, 1] * dx3 + + assert extract_unique_domain_dag(form_scalar) == mesh1 + assert extract_unique_domain_dag(form_vector) == mesh2 + assert extract_unique_domain_dag(form_tensor) == mesh3 + + div_expr = div(u_vector) * f_scalar # Cross-mesh, should fail + with pytest.raises(ValueError, match="Cannot extract unique domain"): + extract_unique_domain_dag(div_expr) + + +def test_extract_unique_domain_repeated_meshes(): + """Test edge cases for domain extraction.""" + cell = triangle + mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=500) + mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=501) + + # MeshSequence with repeated meshes + domain_repeated = MeshSequence([mesh1, mesh2, mesh1]) # mesh1 appears twice + + scalar_elem = LagrangeElement(cell, 1, shape=()) + mixed_elem = MixedElement([scalar_elem, scalar_elem, scalar_elem]) + V = FunctionSpace(domain_repeated, mixed_elem) + u = TrialFunction(V) + + u1, u2, u3 = split(u) + + # Components 0 and 2 should map to mesh1, component 1 to mesh2 + assert extract_unique_domain_dag(u1) == mesh1 # index 0 -> mesh1 + assert extract_unique_domain_dag(u2) == mesh2 # index 1 -> mesh2 + assert extract_unique_domain_dag(u3) == mesh1 # index 2 -> mesh1 (repeated) + + # Expressions combining components on same underlying mesh should work + expr_same = u1 + u3 # Both on mesh1 + assert extract_unique_domain_dag(expr_same) == mesh1 + + # Expressions combining components on different meshes should fail + with pytest.raises(ValueError, match="Cannot extract unique domain"): + extract_unique_domain_dag(u1 + u2) # mesh1 + mesh2 \ No newline at end of file From 062a7791ace37386f14489a292ab17c33ec7d325 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 15:47:54 +0100 Subject: [PATCH 08/35] change `split` to get domain from function space --- ufl/split_functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ufl/split_functions.py b/ufl/split_functions.py index 69be05e92..ecf767100 100644 --- a/ufl/split_functions.py +++ b/ufl/split_functions.py @@ -60,7 +60,7 @@ def split(v): "Don't know how to split tensor valued mixed functions without flattened index space." ) - domain = extract_unique_domain(v, expand_mesh_sequence=False) + domain = v.ufl_function_space().ufl_domain() # Compute value size and set default range end value_size = v.ufl_function_space().value_size From e37249dcb967a12665cfcead8060ec96557613c4 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 15:49:00 +0100 Subject: [PATCH 09/35] add interpolate test --- test/test_unique_domain_extractor.py | 34 +++++++++++++++++++++++++++- 1 file changed, 33 insertions(+), 1 deletion(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index 384a47604..3647fca7e 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -8,6 +8,8 @@ MeshSequence, SpatialCoordinate, Measure, + Argument, + Interpolate, TrialFunction, split, triangle, @@ -261,4 +263,34 @@ def test_extract_unique_domain_repeated_meshes(): # Expressions combining components on different meshes should fail with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(u1 + u2) # mesh1 + mesh2 \ No newline at end of file + extract_unique_domain_dag(u1 + u2) # mesh1 + mesh2 + + +def test_extract_unique_domain_interpolate(): + cell = triangle + mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) + mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) + mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) + domain = MeshSequence([mesh1, mesh2, mesh3]) + scalar_elem = LagrangeElement(cell, 1) # 1 component + vector_elem = LagrangeElement(cell, 1, (2,)) # 2 components + tensor_elem = LagrangeElement(cell, 1, (2, 2)) # 4 components + mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) + V = FunctionSpace(domain, mixed_elem) + + u = TrialFunction(V) + f = Coefficient(V) + + u1, u2, u3 = split(u) + f1, f2, f3 = split(f) + + # Interpolate a function on the mixed space + x1, y1 = SpatialCoordinate(mesh1) + x2, y2 = SpatialCoordinate(mesh2) + x3, y3 = SpatialCoordinate(mesh3) + interp_expr = x1 + cos(u1) * y1 + coarg = Argument(V.dual(), 0) + vstar = split(coarg) + Iu1 = Interpolate(interp_expr, vstar[0]) + expr = Iu1 + f1 * y1 + assert extract_unique_domain_dag(expr) == mesh1 \ No newline at end of file From 85de70dcce2ac94bd05e98da4a2ad18fb056f7d8 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 3 Jul 2025 16:24:11 +0100 Subject: [PATCH 10/35] add test --- test/test_unique_domain_extractor.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index 3647fca7e..f860edb12 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -10,6 +10,7 @@ Measure, Argument, Interpolate, + Constant, TrialFunction, split, triangle, @@ -48,6 +49,9 @@ def test_extract_unique_domain(): expr = u1 + x1 * cos(x1) assert extract_unique_domain_dag(expr) == mesh1 + expr2 = u1 * Constant(mesh1) + x1 + assert extract_unique_domain_dag(expr2) == mesh1 + x2, y2 = SpatialCoordinate(mesh2) with pytest.raises(ValueError) as e_info: _ = extract_unique_domain_dag(u1 + u2) From df1d9116b43bbdfc0367eef245cdba969afb31f2 Mon Sep 17 00:00:00 2001 From: Leo Date: Fri, 4 Jul 2025 12:27:42 +0100 Subject: [PATCH 11/35] add case for `Interpolate` --- ufl/domain.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index b4225e1f6..38cb2dc05 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -26,6 +26,7 @@ from ufl.core.operator import Operator from ufl.core.terminal import Terminal from ufl.indexed import Indexed +from ufl.core.interpolate import Interpolate # Export list for ufl.classes __all_classes__ = ["AbstractDomain", "Mesh", "MeshView"] @@ -334,8 +335,8 @@ def as_domain(domain): (domain,) = set(domain.meshes) return domain try: - return extract_unique_domain(domain) - except AttributeError: + return extract_unique_domain_dag(domain) + except (AttributeError, TypeError): domain = domain.ufl_domain() (domain,) = set(domain.meshes) return domain @@ -557,6 +558,13 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: raise ValueError( f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" ) + + @process.register(Interpolate) + @DAGTraverser.postorder + def _(self, o: Expr, *operand_results) -> AbstractDomain: + """Process Interpolate.""" + fs = o.ufl_function_space() + return fs.ufl_domain() def extract_unique_domain_dag(expr: Union[Expr, Form]) -> AbstractDomain: From 957adaa1f96641e725ec8e16cc2321908f255289 Mon Sep 17 00:00:00 2001 From: Leo Date: Fri, 4 Jul 2025 12:28:02 +0100 Subject: [PATCH 12/35] cyclic import goose chase --- ufl/argument.py | 2 +- ufl/constant.py | 2 +- ufl/form.py | 4 ++-- ufl/functionspace.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/ufl/argument.py b/ufl/argument.py index 043f23dc8..01c9d1445 100644 --- a/ufl/argument.py +++ b/ufl/argument.py @@ -21,7 +21,6 @@ from ufl.duals import is_dual, is_primal from ufl.form import BaseForm from ufl.functionspace import AbstractFunctionSpace, MixedFunctionSpace -from ufl.split_functions import split # Export list for ufl.classes (TODO: not actually classes: drop? these are in ufl.*) __all_classes__ = ["TestFunction", "TrialFunction", "TestFunctions", "TrialFunctions"] @@ -285,6 +284,7 @@ def Arguments(function_space, number): Returns a tuple with the function components corresponding to the subelements. """ + from ufl.split_functions import split if isinstance(function_space, MixedFunctionSpace): return [ Argument(function_space.ufl_sub_space(i), number, i) diff --git a/ufl/constant.py b/ufl/constant.py index dfc65aa06..14959993e 100644 --- a/ufl/constant.py +++ b/ufl/constant.py @@ -8,7 +8,6 @@ from ufl.core.terminal import Terminal from ufl.core.ufl_type import ufl_type -from ufl.domain import as_domain from ufl.utils.counted import Counted @@ -20,6 +19,7 @@ class Constant(Terminal, Counted): def __init__(self, domain, shape=(), count=None): """Initalise.""" + from ufl.domain import as_domain Terminal.__init__(self) Counted.__init__(self, count, Constant) diff --git a/ufl/form.py b/ufl/form.py index b7ff23532..563a98f93 100644 --- a/ufl/form.py +++ b/ufl/form.py @@ -22,7 +22,6 @@ from ufl.core.expr import Expr, ufl_err_str from ufl.core.terminal import FormArgument from ufl.core.ufl_type import UFLType, ufl_type -from ufl.domain import extract_unique_domain, sort_domains from ufl.equation import Equation from ufl.integral import Integral from ufl.utils.counted import Counted @@ -40,6 +39,7 @@ def _sorted_integrals(integrals): Sort integrals by domain id, integral type, subdomain id for a more stable signature computation. """ + from ufl.domain import sort_domains # Group integrals in multilevel dict by keys # [domain][integral_type][subdomain_id] integrals_dict = defaultdict(lambda: defaultdict(lambda: defaultdict(list))) @@ -594,7 +594,7 @@ def __repr__(self): def _analyze_domains(self): """Analyze domains.""" - from ufl.domain import join_domains, sort_domains + from ufl.domain import join_domains, sort_domains, extract_unique_domain # Collect integration domains. self._integration_domains = sort_domains( diff --git a/ufl/functionspace.py b/ufl/functionspace.py index 95757046e..58616ee55 100644 --- a/ufl/functionspace.py +++ b/ufl/functionspace.py @@ -12,7 +12,6 @@ import numpy as np from ufl.core.ufl_type import UFLObject -from ufl.domain import join_domains from ufl.duals import is_dual, is_primal from ufl.utils.sequences import product @@ -311,6 +310,7 @@ def ufl_element(self): def ufl_domains(self): """Return ufl domains.""" + from ufl.domain import join_domains domainlist = [] for s in self._ufl_function_spaces: domainlist.extend(s.ufl_domains()) From 838ae0b6726616cfdcb60cf3d1339fe96d3d5e1d Mon Sep 17 00:00:00 2001 From: Leo Date: Fri, 4 Jul 2025 13:36:30 +0100 Subject: [PATCH 13/35] notes --- ufl/domain.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ufl/domain.py b/ufl/domain.py index 38cb2dc05..7bbc75e14 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -579,6 +579,9 @@ def extract_unique_domain_dag(expr: Union[Expr, Form]) -> AbstractDomain: Returns: AbstractDomain: The unique domain extracted from the expression. """ + # TODO: make this work for BaseForm + # Action - Domain of 0th Argument in result + # Leave AssembledMatrix and Adjoint for now from ufl.form import Form from ufl.core.expr import Expr From 7dffe548679474a55e2e39d82ccc5a50787a6403 Mon Sep 17 00:00:00 2001 From: Leo Date: Tue, 8 Jul 2025 16:32:09 +0100 Subject: [PATCH 14/35] fix imports --- ufl/constant.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ufl/constant.py b/ufl/constant.py index 14959993e..384762550 100644 --- a/ufl/constant.py +++ b/ufl/constant.py @@ -9,6 +9,7 @@ from ufl.core.terminal import Terminal from ufl.core.ufl_type import ufl_type from ufl.utils.counted import Counted +from ufl.domain import as_domain @ufl_type() @@ -19,7 +20,6 @@ class Constant(Terminal, Counted): def __init__(self, domain, shape=(), count=None): """Initalise.""" - from ufl.domain import as_domain Terminal.__init__(self) Counted.__init__(self, count, Constant) From 2de0be5a3943ae0b878f66b8009f4c2f1d894165 Mon Sep 17 00:00:00 2001 From: Leo Date: Wed, 9 Jul 2025 11:34:13 +0100 Subject: [PATCH 15/35] fixes --- ufl/domain.py | 11 +++++++---- ufl/form.py | 3 ++- 2 files changed, 9 insertions(+), 5 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index 7bbc75e14..31eaad67a 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -526,13 +526,16 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: def _(self, o: Expr) -> AbstractDomain: from ufl.coefficient import Coefficient from ufl.argument import Argument + from ufl.geometry import GeometricQuantity + from ufl.constant import Constant if isinstance(o, (Coefficient, Argument)): fs = o.ufl_function_space() return fs.ufl_domain() - try: - return o.ufl_domain() - except AttributeError: - # If the terminal does not have a domain, return None + elif isinstance(o, GeometricQuantity): + return o._domain + elif isinstance(o, Constant): + return o._ufl_domain + else: return None @process.register(Operator) diff --git a/ufl/form.py b/ufl/form.py index 563a98f93..cba128cf9 100644 --- a/ufl/form.py +++ b/ufl/form.py @@ -17,7 +17,7 @@ from itertools import chain from ufl.checks import is_scalar_constant_expression -from ufl.constant import Constant + from ufl.constantvalue import Zero from ufl.core.expr import Expr, ufl_err_str from ufl.core.terminal import FormArgument @@ -399,6 +399,7 @@ def constants(self): def constant_numbering(self): """Return a contiguous numbering of constants in a mapping ``{constant:number}``.""" + from ufl.constant import Constant if self._constant_numbering is None: self._constant_numbering = { expr: num From fa52bf38922434d9d4aa7a1b850fa7b063653f93 Mon Sep 17 00:00:00 2001 From: Leo Date: Tue, 15 Jul 2025 14:03:14 +0100 Subject: [PATCH 16/35] update `measure.py` --- ufl/measure.py | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/ufl/measure.py b/ufl/measure.py index cb9493c64..cba1b3f15 100644 --- a/ufl/measure.py +++ b/ufl/measure.py @@ -15,7 +15,7 @@ from ufl.checks import is_true_ufl_scalar from ufl.constantvalue import as_ufl from ufl.core.expr import Expr -from ufl.domain import AbstractDomain, as_domain, extract_domains +from ufl.domain import AbstractDomain, as_domain, extract_domains, extract_unique_domain_dag from ufl.protocols import id_or_none # Export list for ufl.classes @@ -417,15 +417,16 @@ def __rmul__(self, integrand): # integrand domain = self.ufl_domain() if domain is None: - domains = extract_domains(integrand) - if len(domains) == 1: - (domain,) = domains - elif len(domains) == 0: - raise ValueError("This integral is missing an integration domain.") - else: - raise ValueError( - "Multiple domains found, making the choice of integration domain ambiguous." - ) + domain = extract_unique_domain_dag(integrand) + # domains = extract_domains(integrand) + # if len(domains) == 1: + # (domain,) = domains + # elif len(domains) == 0: + # raise ValueError("This integral is missing an integration domain.") + # else: + # raise ValueError( + # "Multiple domains found, making the choice of integration domain ambiguous." + # ) # Otherwise create and return a one-integral form integral = Integral( From 7befed86889bead75d3f3dbe7f173a67403d4e23 Mon Sep 17 00:00:00 2001 From: Leo Date: Tue, 15 Jul 2025 14:05:35 +0100 Subject: [PATCH 17/35] rename to `unique_domain_extractor` --- test/test_unique_domain_extractor.py | 100 +++++++++++++-------------- ufl/domain.py | 26 +------ ufl/measure.py | 13 +--- 3 files changed, 55 insertions(+), 84 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index f860edb12..b2c8195a7 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -19,7 +19,7 @@ div ) -from ufl.domain import extract_unique_domain_dag +from ufl.domain import extract_unique_domain from ufl.pullback import contravariant_piola, identity_pullback from ufl.sobolevspace import L2, HDiv @@ -38,24 +38,24 @@ def test_extract_unique_domain(): u = TrialFunction(V) u1, u2, u3 = split(u) for i, u_i in enumerate((u1, u2, u3)): - assert extract_unique_domain_dag(u_i) == domain[i] + assert extract_unique_domain(u_i) == domain[i] f = Coefficient(V) f1, f2, f3 = split(f) for i, f_i in enumerate((f1, f2, f3)): - assert extract_unique_domain_dag(f_i) == domain[i] + assert extract_unique_domain(f_i) == domain[i] x1, y1 = SpatialCoordinate(mesh1) expr = u1 + x1 * cos(x1) - assert extract_unique_domain_dag(expr) == mesh1 + assert extract_unique_domain(expr) == mesh1 expr2 = u1 * Constant(mesh1) + x1 - assert extract_unique_domain_dag(expr2) == mesh1 + assert extract_unique_domain(expr2) == mesh1 x2, y2 = SpatialCoordinate(mesh2) with pytest.raises(ValueError) as e_info: - _ = extract_unique_domain_dag(u1 + u2) - _ = extract_unique_domain_dag(u1 + u2 + x2 * cos(x2 * u1)) + _ = extract_unique_domain(u1 + u2) + _ = extract_unique_domain(u1 + u2 + x2 * cos(x2 * u1)) def test_extract_unique_domain_form(): @@ -80,7 +80,7 @@ def test_extract_unique_domain_form(): form1 = inner(u1, f1) * dx1 - assert extract_unique_domain_dag(form1) == mesh1 + assert extract_unique_domain(form1) == mesh1 def test_extract_unique_domain_single_mesh(): @@ -95,8 +95,8 @@ def test_extract_unique_domain_single_mesh(): f_scalar = Coefficient(V_scalar) # Test that single mesh functions return the correct domain - assert extract_unique_domain_dag(u_scalar) == mesh - assert extract_unique_domain_dag(f_scalar) == mesh + assert extract_unique_domain(u_scalar) == mesh + assert extract_unique_domain(f_scalar) == mesh # Test vector elements P1_vec = LagrangeElement(cell, 1, (2,)) @@ -104,14 +104,14 @@ def test_extract_unique_domain_single_mesh(): u_vector = TrialFunction(V_vector) f_vector = Coefficient(V_vector) - assert extract_unique_domain_dag(u_vector) == mesh - assert extract_unique_domain_dag(f_vector) == mesh + assert extract_unique_domain(u_vector) == mesh + assert extract_unique_domain(f_vector) == mesh # Test indexing into vector elements - assert extract_unique_domain_dag(u_vector[0]) == mesh - assert extract_unique_domain_dag(u_vector[1]) == mesh - assert extract_unique_domain_dag(f_vector[0]) == mesh - assert extract_unique_domain_dag(f_vector[1]) == mesh + assert extract_unique_domain(u_vector[0]) == mesh + assert extract_unique_domain(u_vector[1]) == mesh + assert extract_unique_domain(f_vector[0]) == mesh + assert extract_unique_domain(f_vector[1]) == mesh # Test tensor elements P1_tensor = LagrangeElement(cell, 1, (2, 2)) @@ -119,11 +119,11 @@ def test_extract_unique_domain_single_mesh(): u_tensor = TrialFunction(V_tensor) f_tensor = Coefficient(V_tensor) - assert extract_unique_domain_dag(u_tensor) == mesh - assert extract_unique_domain_dag(f_tensor) == mesh - assert extract_unique_domain_dag(u_tensor[0, 0]) == mesh - assert extract_unique_domain_dag(u_tensor[1, 1]) == mesh - assert extract_unique_domain_dag(f_tensor[0, 1]) == mesh + assert extract_unique_domain(u_tensor) == mesh + assert extract_unique_domain(f_tensor) == mesh + assert extract_unique_domain(u_tensor[0, 0]) == mesh + assert extract_unique_domain(u_tensor[1, 1]) == mesh + assert extract_unique_domain(f_tensor[0, 1]) == mesh # Test expressions combining same-domain functions x, y = SpatialCoordinate(mesh) @@ -131,14 +131,14 @@ def test_extract_unique_domain_single_mesh(): expr2 = u_vector[0] + x expr3 = inner(u_vector, f_vector) - assert extract_unique_domain_dag(expr1) == mesh - assert extract_unique_domain_dag(expr2) == mesh - assert extract_unique_domain_dag(expr3) == mesh + assert extract_unique_domain(expr1) == mesh + assert extract_unique_domain(expr2) == mesh + assert extract_unique_domain(expr3) == mesh # Test forms dx = Measure("dx", mesh) form = inner(u_scalar, f_scalar) * dx - assert extract_unique_domain_dag(form) == mesh + assert extract_unique_domain(form) == mesh def test_extract_unique_domain_mixed_scalar_vector_tensor(): @@ -171,19 +171,19 @@ def test_extract_unique_domain_mixed_scalar_vector_tensor(): # Test that each component maps to correct mesh for i, u_i in enumerate((u_scalar, u_vector, u_tensor)): - assert extract_unique_domain_dag(u_i) == domain[i] + assert extract_unique_domain(u_i) == domain[i] for i, f_i in enumerate((f_scalar, f_vector, f_tensor)): - assert extract_unique_domain_dag(f_i) == domain[i] + assert extract_unique_domain(f_i) == domain[i] # Test indexing into vector and tensor components for i in range(2): - assert extract_unique_domain_dag(u_vector[i]) == mesh2 - assert extract_unique_domain_dag(f_vector[i]) == mesh2 + assert extract_unique_domain(u_vector[i]) == mesh2 + assert extract_unique_domain(f_vector[i]) == mesh2 for i in range(2): for j in range(2): - assert extract_unique_domain_dag(u_tensor[i, j]) == mesh3 - assert extract_unique_domain_dag(f_tensor[i, j]) == mesh3 + assert extract_unique_domain(u_tensor[i, j]) == mesh3 + assert extract_unique_domain(f_tensor[i, j]) == mesh3 # Test expressions on same mesh (should work) x1, y1 = SpatialCoordinate(mesh1) @@ -192,35 +192,35 @@ def test_extract_unique_domain_mixed_scalar_vector_tensor(): # Scalar expressions expr_scalar = u_scalar * y1 + f_scalar + x1 - assert extract_unique_domain_dag(expr_scalar) == mesh1 + assert extract_unique_domain(expr_scalar) == mesh1 # Vector expressions expr_vector = inner(u_vector * y2, f_vector) + x2 - assert extract_unique_domain_dag(expr_vector) == mesh2 + assert extract_unique_domain(expr_vector) == mesh2 # Vector component expressions expr_vec_comp = u_vector[0] + f_vector[1] * y2 + x2 - assert extract_unique_domain_dag(expr_vec_comp) == mesh2 + assert extract_unique_domain(expr_vec_comp) == mesh2 # Tensor expressions expr_tensor = y3 * u_tensor[0, 0] + f_tensor[1, 1] + x3 - assert extract_unique_domain_dag(expr_tensor) == mesh3 + assert extract_unique_domain(expr_tensor) == mesh3 # Test expressions mixing different mesh components (should fail) with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(u_scalar + u_vector[0]) + extract_unique_domain(u_scalar + u_vector[0]) with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(u_vector[0] + u_tensor[0, 0]) + extract_unique_domain(u_vector[0] + u_tensor[0, 0]) with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(f_scalar + f_tensor[1, 1]) + extract_unique_domain(f_scalar + f_tensor[1, 1]) with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(u_scalar + x2) + extract_unique_domain(u_scalar + x2) with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(u_vector[0] + x3) + extract_unique_domain(u_vector[0] + x3) # Test forms on individual meshes dx1 = Measure("dx", mesh1) @@ -231,13 +231,13 @@ def test_extract_unique_domain_mixed_scalar_vector_tensor(): form_vector = inner(u_vector, f_vector) * dx2 form_tensor = u_tensor[0, 0] * f_tensor[1, 1] * dx3 - assert extract_unique_domain_dag(form_scalar) == mesh1 - assert extract_unique_domain_dag(form_vector) == mesh2 - assert extract_unique_domain_dag(form_tensor) == mesh3 + assert extract_unique_domain(form_scalar) == mesh1 + assert extract_unique_domain(form_vector) == mesh2 + assert extract_unique_domain(form_tensor) == mesh3 div_expr = div(u_vector) * f_scalar # Cross-mesh, should fail with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(div_expr) + extract_unique_domain(div_expr) def test_extract_unique_domain_repeated_meshes(): @@ -257,17 +257,17 @@ def test_extract_unique_domain_repeated_meshes(): u1, u2, u3 = split(u) # Components 0 and 2 should map to mesh1, component 1 to mesh2 - assert extract_unique_domain_dag(u1) == mesh1 # index 0 -> mesh1 - assert extract_unique_domain_dag(u2) == mesh2 # index 1 -> mesh2 - assert extract_unique_domain_dag(u3) == mesh1 # index 2 -> mesh1 (repeated) + assert extract_unique_domain(u1) == mesh1 # index 0 -> mesh1 + assert extract_unique_domain(u2) == mesh2 # index 1 -> mesh2 + assert extract_unique_domain(u3) == mesh1 # index 2 -> mesh1 (repeated) # Expressions combining components on same underlying mesh should work expr_same = u1 + u3 # Both on mesh1 - assert extract_unique_domain_dag(expr_same) == mesh1 + assert extract_unique_domain(expr_same) == mesh1 # Expressions combining components on different meshes should fail with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain_dag(u1 + u2) # mesh1 + mesh2 + extract_unique_domain(u1 + u2) # mesh1 + mesh2 def test_extract_unique_domain_interpolate(): @@ -297,4 +297,4 @@ def test_extract_unique_domain_interpolate(): vstar = split(coarg) Iu1 = Interpolate(interp_expr, vstar[0]) expr = Iu1 + f1 * y1 - assert extract_unique_domain_dag(expr) == mesh1 \ No newline at end of file + assert extract_unique_domain(expr) == mesh1 \ No newline at end of file diff --git a/ufl/domain.py b/ufl/domain.py index 31eaad67a..c27ef0e89 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -335,7 +335,7 @@ def as_domain(domain): (domain,) = set(domain.meshes) return domain try: - return extract_unique_domain_dag(domain) + return extract_unique_domain(domain) except (AttributeError, TypeError): domain = domain.ufl_domain() (domain,) = set(domain.meshes) @@ -416,26 +416,6 @@ def extract_domains(expr: Expr | Form, expand_mesh_sequence: bool = True): return sort_domains(join_domains(domainlist, expand_mesh_sequence=expand_mesh_sequence)) -def extract_unique_domain(expr, expand_mesh_sequence: bool = True): - """Return the single unique domain expression is defined on or throw an error. - - Args: - expr: Expr or Form. - expand_mesh_sequence: If True, MeshSequence components are expanded. - - Returns: - domain. - - """ - domains = extract_domains(expr, expand_mesh_sequence=expand_mesh_sequence) - if len(domains) == 1: - return domains[0] - elif domains: - raise ValueError("Found multiple domains, cannot return just one.") - else: - return None - - def find_geometric_dimension(expr): """Find the geometric dimension of an expression.""" gdims = set() @@ -570,7 +550,7 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: return fs.ufl_domain() -def extract_unique_domain_dag(expr: Union[Expr, Form]) -> AbstractDomain: +def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: """Extract the single unique domain from an expression. This works for expressions containing Indexed Arguments and Coefficients from @@ -592,7 +572,7 @@ def extract_unique_domain_dag(expr: Union[Expr, Form]) -> AbstractDomain: # For forms, we extract domains from integrals domains = set() for integral in expr.integrals(): - domain = extract_unique_domain_dag(integral.integrand()) + domain = extract_unique_domain(integral.integrand()) if domain is not None: domains.add(domain) diff --git a/ufl/measure.py b/ufl/measure.py index cba1b3f15..27415b85e 100644 --- a/ufl/measure.py +++ b/ufl/measure.py @@ -15,7 +15,7 @@ from ufl.checks import is_true_ufl_scalar from ufl.constantvalue import as_ufl from ufl.core.expr import Expr -from ufl.domain import AbstractDomain, as_domain, extract_domains, extract_unique_domain_dag +from ufl.domain import AbstractDomain, as_domain, extract_unique_domain from ufl.protocols import id_or_none # Export list for ufl.classes @@ -417,16 +417,7 @@ def __rmul__(self, integrand): # integrand domain = self.ufl_domain() if domain is None: - domain = extract_unique_domain_dag(integrand) - # domains = extract_domains(integrand) - # if len(domains) == 1: - # (domain,) = domains - # elif len(domains) == 0: - # raise ValueError("This integral is missing an integration domain.") - # else: - # raise ValueError( - # "Multiple domains found, making the choice of integration domain ambiguous." - # ) + domain = extract_unique_domain(integrand) # Otherwise create and return a one-integral form integral = Integral( From b3770dfc2b9301b0d634658319d2bf9d84c4eaad Mon Sep 17 00:00:00 2001 From: Leo Date: Tue, 15 Jul 2025 14:11:09 +0100 Subject: [PATCH 18/35] remove relevant uses of `expand_mesh_sequence` --- ufl/algorithms/apply_derivatives.py | 4 ++-- ufl/algorithms/compute_form_data.py | 4 ++-- ufl/differentiation.py | 4 ++-- ufl/form.py | 2 +- ufl/pullback.py | 4 ++-- 5 files changed, 9 insertions(+), 9 deletions(-) diff --git a/ufl/algorithms/apply_derivatives.py b/ufl/algorithms/apply_derivatives.py index 3c5a07919..f29eae9ac 100644 --- a/ufl/algorithms/apply_derivatives.py +++ b/ufl/algorithms/apply_derivatives.py @@ -834,7 +834,7 @@ def _(self, o: ReferenceValue) -> Expr: f = o.ufl_operands[0] if not f._ufl_is_terminal_: raise ValueError("ReferenceValue can only wrap a terminal") - domain = extract_unique_domain(f, expand_mesh_sequence=False) + domain = extract_unique_domain(f) if isinstance(domain, MeshSequence): element = f.ufl_function_space().ufl_element() # type: ignore if element.num_sub_elements != len(domain): @@ -897,7 +897,7 @@ def _(self, o: Expr) -> Expr: ) if not valid_operand: raise ValueError("ReferenceGrad can only wrap a reference frame type!") - domain = extract_unique_domain(f, expand_mesh_sequence=False) + domain = extract_unique_domain(f) if isinstance(domain, MeshSequence): if not f._ufl_is_in_reference_frame_: raise RuntimeError("Expecting a reference frame type") diff --git a/ufl/algorithms/compute_form_data.py b/ufl/algorithms/compute_form_data.py index 09f98af4e..10829eeed 100644 --- a/ufl/algorithms/compute_form_data.py +++ b/ufl/algorithms/compute_form_data.py @@ -188,7 +188,7 @@ def _build_coefficient_replace_map(coefficients, element_mapping=None): # coefficient had a domain, the new one does too. # This should be overhauled with requirement that Expressions # always have a domain. - domain = extract_unique_domain(f, expand_mesh_sequence=False) + domain = extract_unique_domain(f) if domain is not None: new_e = FunctionSpace(domain, new_e) new_f = Coefficient(new_e, count=i) @@ -454,7 +454,7 @@ def compute_form_data( for o in self.reduced_coefficients: if o in coefficients_to_split: c = self.function_replace_map[o] - mesh = extract_unique_domain(c, expand_mesh_sequence=False) + mesh = extract_unique_domain(c) elem = c.ufl_element() coefficient_split[c] = [ Coefficient(FunctionSpace(m, e)) diff --git a/ufl/differentiation.py b/ufl/differentiation.py index fd4a66bb7..76d327f8e 100644 --- a/ufl/differentiation.py +++ b/ufl/differentiation.py @@ -318,7 +318,7 @@ def __new__(cls, f): # Return zero if expression is trivially constant if is_cellwise_constant(f): # TODO: Use max topological dimension if there are multiple topological dimensions. - dim = extract_unique_domain(f, expand_mesh_sequence=False).topological_dimension() + dim = extract_unique_domain(f).topological_dimension() return Zero(f.ufl_shape + (dim,), f.ufl_free_indices, f.ufl_index_dimensions) return CompoundDerivative.__new__(cls) @@ -326,7 +326,7 @@ def __init__(self, f): """Initalise.""" CompoundDerivative.__init__(self, (f,)) # TODO: Use max topological dimension if there are multiple topological dimensions. - self._dim = extract_unique_domain(f, expand_mesh_sequence=False).topological_dimension() + self._dim = extract_unique_domain(f).topological_dimension() def _ufl_expr_reconstruct_(self, op): """Return a new object of the same type with new operands.""" diff --git a/ufl/form.py b/ufl/form.py index cba128cf9..3f5ad0c17 100644 --- a/ufl/form.py +++ b/ufl/form.py @@ -606,7 +606,7 @@ def _analyze_domains(self): for o in chain( self.arguments(), self.coefficients(), self.constants(), self.geometric_quantities() ): - domain = extract_unique_domain(o, expand_mesh_sequence=False) + domain = extract_unique_domain(o) domains_in_integrands.update(domain.meshes) domains_in_integrands -= set(self._integration_domains) all_domains = self._integration_domains + sort_domains(join_domains(domains_in_integrands)) diff --git a/ufl/pullback.py b/ufl/pullback.py index 4f4bacdc6..fa2f1362b 100644 --- a/ufl/pullback.py +++ b/ufl/pullback.py @@ -414,7 +414,7 @@ def apply(self, expr, domain=None): g_components = [] offset = 0 # For each unique piece in reference space, apply the appropriate pullback - domain = domain or extract_unique_domain(expr, expand_mesh_sequence=False) + domain = domain or extract_unique_domain(expr) if isinstance(domain, MeshSequence): if len(domain) != self._element.num_sub_elements: raise ValueError(f"""num. component meshes ({len(domain)}) != @@ -504,7 +504,7 @@ def apply(self, expr, domain=None): for subelem in self._element.sub_elements: offsets.append(offsets[-1] + subelem.reference_value_size) # For each unique piece in reference space, apply the appropriate pullback - domain = domain or extract_unique_domain(expr, expand_mesh_sequence=False) + domain = domain or extract_unique_domain(expr) if isinstance(domain, MeshSequence): if len(domain) != self._element.num_sub_elements: raise ValueError(f"""num. component meshes ({len(domain)}) != From f5490bafd9829056be38e08501241a3a03182d6a Mon Sep 17 00:00:00 2001 From: Leo Date: Tue, 15 Jul 2025 14:14:43 +0100 Subject: [PATCH 19/35] tidy up tidy up update comment --- ufl/domain.py | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index c27ef0e89..523b16160 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -433,9 +433,6 @@ def find_geometric_dimension(expr): class UniqueDomainExtractor(DAGTraverser): - # Terminals: Arg or Coeff or constant. constant has no domain. - # Operators: 1) Children match domain, return that. 2) Differing domains, raise error. 3) Different domains, return A or B - # ExternalOperator or Interpolate: Domain is where we are mapping into. These are both BaseFormOperators. def __init__( self, @@ -473,31 +470,26 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: # Get the domain from the first operand (the expression being indexed) expression_domain = operand_results[0] - # If the expression domain is a MeshSequence, extract the specific mesh if isinstance(expression_domain, MeshSequence): index = multiindex[0]._value - # For mixed function spaces with MeshSequence, we need to determine - # which component mesh this index corresponds to - # This requires understanding the mixed element structure element = expression.ufl_element() if hasattr(element, 'sub_elements'): - # Calculate which sub-element (and thus which mesh) this index belongs to + # Need to do this in case we have sub elements which are vector or tensor valued offset = 0 for i, sub_element in enumerate(element.sub_elements): # Get the value size for this sub-element on its corresponding mesh - mesh_for_element = expression_domain.meshes[i] - sub_element_fs = FunctionSpace(mesh_for_element, sub_element) + sub_element_mesh = expression_domain.meshes[i] + sub_element_fs = FunctionSpace(sub_element_mesh, sub_element) sub_element_size = sub_element_fs.value_size if index < offset + sub_element_size: # This index belongs to mesh i - return mesh_for_element + return sub_element_mesh offset += sub_element_size raise ValueError(f"Index {index} out of range for mixed function space") else: - # Simple case: direct indexing into MeshSequence - if 0 <= index < len(expression_domain.meshes): - return expression_domain.meshes[index] + # If no sub elements we just grab the mesh + return expression_domain.meshes[index] # If it's not a MeshSequence, just return the expression domain return expression_domain @@ -537,7 +529,7 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: if all(d == first_domain for d in domains): return first_domain else: - # If domains are not none and differ, raise error + # If domains are not none and differ, then we don't have a unique domain raise ValueError( f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" ) @@ -569,7 +561,6 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: from ufl.core.expr import Expr if isinstance(expr, Form): - # For forms, we extract domains from integrals domains = set() for integral in expr.integrals(): domain = extract_unique_domain(integral.integrand()) @@ -583,8 +574,7 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: else: raise ValueError(f"Form has multiple domains: {domains}") elif isinstance(expr, Expr): - # For expressions, use the DAG traverser extractor = UniqueDomainExtractor() return extractor(expr) else: - raise TypeError(f"Unsupported type for domain extraction: {type(expr)}. Expected Expr or Form.") \ No newline at end of file + raise TypeError(f"Unsupported type for domain extraction: {type(expr).__name__}. Expected Expr or Form.") From 8555dfa8add8b7aa382f8c172145895d2ccc79f6 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 17 Jul 2025 12:16:18 +0100 Subject: [PATCH 20/35] lint lint --- ufl/domain.py | 30 +++++++++++++++--------------- ufl/split_functions.py | 1 - 2 files changed, 15 insertions(+), 16 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index 523b16160..3dd8d3657 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -458,7 +458,7 @@ def process(self, o: Expr) -> Expr: """ return super().process(o) - + @process.register(Indexed) @DAGTraverser.postorder def _(self, o: Expr, *operand_results) -> AbstractDomain: @@ -466,10 +466,10 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: from ufl.functionspace import FunctionSpace # Get the indexed expression and the multiindex expression, multiindex = o.ufl_operands - + # Get the domain from the first operand (the expression being indexed) expression_domain = operand_results[0] - + if isinstance(expression_domain, MeshSequence): index = multiindex[0]._value element = expression.ufl_element() @@ -481,7 +481,7 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: sub_element_mesh = expression_domain.meshes[i] sub_element_fs = FunctionSpace(sub_element_mesh, sub_element) sub_element_size = sub_element_fs.value_size - + if index < offset + sub_element_size: # This index belongs to mesh i return sub_element_mesh @@ -490,7 +490,7 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: else: # If no sub elements we just grab the mesh return expression_domain.meshes[index] - + # If it's not a MeshSequence, just return the expression domain return expression_domain @@ -509,14 +509,14 @@ def _(self, o: Expr) -> AbstractDomain: return o._ufl_domain else: return None - + @process.register(Operator) @DAGTraverser.postorder def _(self, o: Expr, *operand_results) -> AbstractDomain: """Process Operator.""" # Filter out None results (from operands that don't have domains) domains = [d for d in operand_results if d is not None] - + if not domains: # No operands have domains return None @@ -533,24 +533,24 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: raise ValueError( f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" ) - + @process.register(Interpolate) @DAGTraverser.postorder def _(self, o: Expr, *operand_results) -> AbstractDomain: """Process Interpolate.""" fs = o.ufl_function_space() - return fs.ufl_domain() + return fs.ufl_domain() def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: """Extract the single unique domain from an expression. - + This works for expressions containing Indexed Arguments and Coefficients from split functions on mixed function spaces. - + Args: expr: Expr or Form to extract domain from - + Returns: AbstractDomain: The unique domain extracted from the expression. """ @@ -559,14 +559,14 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: # Leave AssembledMatrix and Adjoint for now from ufl.form import Form from ufl.core.expr import Expr - + if isinstance(expr, Form): domains = set() for integral in expr.integrals(): domain = extract_unique_domain(integral.integrand()) if domain is not None: domains.add(domain) - + if len(domains) == 0: return None elif len(domains) == 1: @@ -577,4 +577,4 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: extractor = UniqueDomainExtractor() return extractor(expr) else: - raise TypeError(f"Unsupported type for domain extraction: {type(expr).__name__}. Expected Expr or Form.") + raise TypeError(f"Expected an Expr or Form, not a {type(expr).__name__}.") diff --git a/ufl/split_functions.py b/ufl/split_functions.py index ecf767100..5ece9d262 100644 --- a/ufl/split_functions.py +++ b/ufl/split_functions.py @@ -7,7 +7,6 @@ # # Modified by Anders Logg, 2008 -from ufl.domain import extract_unique_domain from ufl.functionspace import FunctionSpace from ufl.indexed import Indexed from ufl.permutation import compute_indices From 5cb1c574703c0d5802d83037c853933dc5d91b59 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 17 Jul 2025 13:45:04 +0100 Subject: [PATCH 21/35] base form --- ufl/domain.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/ufl/domain.py b/ufl/domain.py index 3dd8d3657..c7d7aed14 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -557,7 +557,7 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: # TODO: make this work for BaseForm # Action - Domain of 0th Argument in result # Leave AssembledMatrix and Adjoint for now - from ufl.form import Form + from ufl.form import Form, BaseForm from ufl.core.expr import Expr if isinstance(expr, Form): @@ -573,6 +573,8 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: return domains.pop() else: raise ValueError(f"Form has multiple domains: {domains}") + elif isinstance(expr, BaseForm): + return extract_unique_domain(expr.arguments()[0]) elif isinstance(expr, Expr): extractor = UniqueDomainExtractor() return extractor(expr) From 1c77207dee77319ecf28aa98e87488c3f6c22e43 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 17 Jul 2025 13:45:29 +0100 Subject: [PATCH 22/35] update tests fixes --- test/test_unique_domain_extractor.py | 192 +++++++++++++++------------ 1 file changed, 104 insertions(+), 88 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index b2c8195a7..042f0df67 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -12,11 +12,19 @@ Interpolate, Constant, TrialFunction, + TrialFunctions, + TestFunction, + TestFunctions, split, triangle, cos, inner, - div + div, + FacetNormal, + grad, + action, + MixedFunctionSpace, + Matrix ) from ufl.domain import extract_unique_domain @@ -72,16 +80,27 @@ def test_extract_unique_domain_form(): u = TrialFunction(V) u1, u2, u3 = split(u) + v = TestFunction(V) + v1, v2, v3 = split(v) f = Coefficient(V) f1, f2, f3 = split(f) + n = FacetNormal(mesh1) dx1 = Measure("dx", mesh1) + ds1 = Measure("ds", mesh1) + dx2 = Measure("dx", mesh2) - form1 = inner(u1, f1) * dx1 - + form1 = inner(grad(u1), grad(v1)) * dx1 - inner(grad(u1), n) * v1 * ds1 assert extract_unique_domain(form1) == mesh1 + form2 = inner(u1, f1) * dx1 + assert extract_unique_domain(form2) == mesh1 + + form3 = inner(u1, v1) * dx1 + inner(u2, v2) * dx2 + with pytest.raises(ValueError): + extract_unique_domain(form3) + def test_extract_unique_domain_single_mesh(): """Test domain extraction for standard function spaces on a single mesh.""" @@ -94,11 +113,9 @@ def test_extract_unique_domain_single_mesh(): u_scalar = TrialFunction(V_scalar) f_scalar = Coefficient(V_scalar) - # Test that single mesh functions return the correct domain assert extract_unique_domain(u_scalar) == mesh assert extract_unique_domain(f_scalar) == mesh - # Test vector elements P1_vec = LagrangeElement(cell, 1, (2,)) V_vector = FunctionSpace(mesh, P1_vec) u_vector = TrialFunction(V_vector) @@ -107,13 +124,11 @@ def test_extract_unique_domain_single_mesh(): assert extract_unique_domain(u_vector) == mesh assert extract_unique_domain(f_vector) == mesh - # Test indexing into vector elements assert extract_unique_domain(u_vector[0]) == mesh assert extract_unique_domain(u_vector[1]) == mesh assert extract_unique_domain(f_vector[0]) == mesh assert extract_unique_domain(f_vector[1]) == mesh - # Test tensor elements P1_tensor = LagrangeElement(cell, 1, (2, 2)) V_tensor = FunctionSpace(mesh, P1_tensor) u_tensor = TrialFunction(V_tensor) @@ -125,7 +140,6 @@ def test_extract_unique_domain_single_mesh(): assert extract_unique_domain(u_tensor[1, 1]) == mesh assert extract_unique_domain(f_tensor[0, 1]) == mesh - # Test expressions combining same-domain functions x, y = SpatialCoordinate(mesh) expr1 = u_scalar + f_scalar expr2 = u_vector[0] + x @@ -142,40 +156,31 @@ def test_extract_unique_domain_single_mesh(): def test_extract_unique_domain_mixed_scalar_vector_tensor(): - """Test domain extraction for mixed function spaces with scalar, vector, and tensor elements.""" + """Test domain extraction for mixed function spaces + with scalar, vector, and tensor elements.""" cell = triangle mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) domain = MeshSequence([mesh1, mesh2, mesh3]) - - # Create mixed element with different types: - # - scalar on mesh1 - # - vector on mesh2 - # - tensor on mesh3 - scalar_elem = LagrangeElement(cell, 1) # 1 component - vector_elem = LagrangeElement(cell, 1, (2,)) # 2 components - tensor_elem = LagrangeElement(cell, 1, (2, 2)) # 4 components + + scalar_elem = LagrangeElement(cell, 1) + vector_elem = LagrangeElement(cell, 1, (2,)) + tensor_elem = LagrangeElement(cell, 1, (2, 2)) mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) - + V = FunctionSpace(domain, mixed_elem) u = TrialFunction(V) f = Coefficient(V) - - # For mixed element [scalar, vector(2), tensor(2,2)], split gives: - # - u_components[0]: scalar (mesh1) - index 0 - # - u_components[1]: vector (mesh2) - indices 1-2 - # - u_components[2]: tensor (mesh3) - indices 3-6 + u_scalar, u_vector, u_tensor = split(u) f_scalar, f_vector, f_tensor = split(f) - - # Test that each component maps to correct mesh + for i, u_i in enumerate((u_scalar, u_vector, u_tensor)): assert extract_unique_domain(u_i) == domain[i] for i, f_i in enumerate((f_scalar, f_vector, f_tensor)): assert extract_unique_domain(f_i) == domain[i] - # Test indexing into vector and tensor components for i in range(2): assert extract_unique_domain(u_vector[i]) == mesh2 assert extract_unique_domain(f_vector[i]) == mesh2 @@ -185,116 +190,127 @@ def test_extract_unique_domain_mixed_scalar_vector_tensor(): assert extract_unique_domain(u_tensor[i, j]) == mesh3 assert extract_unique_domain(f_tensor[i, j]) == mesh3 - # Test expressions on same mesh (should work) x1, y1 = SpatialCoordinate(mesh1) x2, y2 = SpatialCoordinate(mesh2) x3, y3 = SpatialCoordinate(mesh3) - - # Scalar expressions + expr_scalar = u_scalar * y1 + f_scalar + x1 assert extract_unique_domain(expr_scalar) == mesh1 - - # Vector expressions + expr_vector = inner(u_vector * y2, f_vector) + x2 assert extract_unique_domain(expr_vector) == mesh2 - - # Vector component expressions + expr_vec_comp = u_vector[0] + f_vector[1] * y2 + x2 assert extract_unique_domain(expr_vec_comp) == mesh2 - - # Tensor expressions + expr_tensor = y3 * u_tensor[0, 0] + f_tensor[1, 1] + x3 assert extract_unique_domain(expr_tensor) == mesh3 - - # Test expressions mixing different mesh components (should fail) - with pytest.raises(ValueError, match="Cannot extract unique domain"): + + with pytest.raises(ValueError): extract_unique_domain(u_scalar + u_vector[0]) - - with pytest.raises(ValueError, match="Cannot extract unique domain"): + + with pytest.raises(ValueError): extract_unique_domain(u_vector[0] + u_tensor[0, 0]) - - with pytest.raises(ValueError, match="Cannot extract unique domain"): + + with pytest.raises(ValueError): extract_unique_domain(f_scalar + f_tensor[1, 1]) - - with pytest.raises(ValueError, match="Cannot extract unique domain"): + + with pytest.raises(ValueError): extract_unique_domain(u_scalar + x2) - - with pytest.raises(ValueError, match="Cannot extract unique domain"): + + with pytest.raises(ValueError): extract_unique_domain(u_vector[0] + x3) - - # Test forms on individual meshes + dx1 = Measure("dx", mesh1) dx2 = Measure("dx", mesh2) dx3 = Measure("dx", mesh3) - + form_scalar = u_scalar * f_scalar * dx1 form_vector = inner(u_vector, f_vector) * dx2 form_tensor = u_tensor[0, 0] * f_tensor[1, 1] * dx3 - + assert extract_unique_domain(form_scalar) == mesh1 assert extract_unique_domain(form_vector) == mesh2 assert extract_unique_domain(form_tensor) == mesh3 - - div_expr = div(u_vector) * f_scalar # Cross-mesh, should fail - with pytest.raises(ValueError, match="Cannot extract unique domain"): + + div_expr = div(u_vector) * f_scalar + with pytest.raises(ValueError): extract_unique_domain(div_expr) def test_extract_unique_domain_repeated_meshes(): - """Test edge cases for domain extraction.""" cell = triangle mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=500) mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=501) - + # MeshSequence with repeated meshes - domain_repeated = MeshSequence([mesh1, mesh2, mesh1]) # mesh1 appears twice - + domain_repeated = MeshSequence([mesh1, mesh2, mesh1]) + scalar_elem = LagrangeElement(cell, 1, shape=()) mixed_elem = MixedElement([scalar_elem, scalar_elem, scalar_elem]) V = FunctionSpace(domain_repeated, mixed_elem) u = TrialFunction(V) - + u1, u2, u3 = split(u) + + assert extract_unique_domain(u1) == mesh1 + assert extract_unique_domain(u2) == mesh2 + assert extract_unique_domain(u3) == mesh1 + + expr_same = u1 + u3 + assert extract_unique_domain(expr_same) == mesh1 + + with pytest.raises(ValueError): + extract_unique_domain(u1 + u2) + + +def test_extract_unique_domain_action(): + cell = triangle + mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) + mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) + mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) + domain = MeshSequence([mesh1, mesh2, mesh3]) + scalar_elem = LagrangeElement(cell, 1) + vector_elem = LagrangeElement(cell, 1, (2,)) + tensor_elem = LagrangeElement(cell, 1, (2, 2)) + mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) + + V1 = FunctionSpace(mesh1, scalar_elem) + V2 = FunctionSpace(mesh2, scalar_elem) - # Components 0 and 2 should map to mesh1, component 1 to mesh2 - assert extract_unique_domain(u1) == mesh1 # index 0 -> mesh1 - assert extract_unique_domain(u2) == mesh2 # index 1 -> mesh2 - assert extract_unique_domain(u3) == mesh1 # index 2 -> mesh1 (repeated) + A = Matrix(V1, V2) + u = Coefficient(V2) - # Expressions combining components on same underlying mesh should work - expr_same = u1 + u3 # Both on mesh1 - assert extract_unique_domain(expr_same) == mesh1 + action_Au = action(A, u) + assert extract_unique_domain(action_Au) == mesh1 - # Expressions combining components on different meshes should fail - with pytest.raises(ValueError, match="Cannot extract unique domain"): - extract_unique_domain(u1 + u2) # mesh1 + mesh2 + V_mixed1 = FunctionSpace(domain, mixed_elem) + V_mixed2 = FunctionSpace(domain, mixed_elem) + A_mixed = Matrix(V_mixed1, V_mixed2) + u_mixed = Coefficient(V_mixed2) + + action_mixed = action(A_mixed, u_mixed) + assert extract_unique_domain(action_mixed) == domain + def test_extract_unique_domain_interpolate(): cell = triangle mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) - domain = MeshSequence([mesh1, mesh2, mesh3]) - scalar_elem = LagrangeElement(cell, 1) # 1 component - vector_elem = LagrangeElement(cell, 1, (2,)) # 2 components - tensor_elem = LagrangeElement(cell, 1, (2, 2)) # 4 components - mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) - V = FunctionSpace(domain, mixed_elem) + scalar_elem = LagrangeElement(cell, 1) + vector_elem = LagrangeElement(cell, 1, (2,)) + tensor_elem = LagrangeElement(cell, 1, (2, 2)) - u = TrialFunction(V) - f = Coefficient(V) + V_1 = FunctionSpace(mesh1, scalar_elem) + V_2 = FunctionSpace(mesh2, vector_elem) + V_3 = FunctionSpace(mesh3, tensor_elem) - u1, u2, u3 = split(u) - f1, f2, f3 = split(f) + # # MixedFunctionSpace = V_3d x V_2d x V_1d + # V = MixedFunctionSpace(V_1, V_2, V_3) - # Interpolate a function on the mixed space - x1, y1 = SpatialCoordinate(mesh1) - x2, y2 = SpatialCoordinate(mesh2) - x3, y3 = SpatialCoordinate(mesh3) - interp_expr = x1 + cos(u1) * y1 - coarg = Argument(V.dual(), 0) - vstar = split(coarg) - Iu1 = Interpolate(interp_expr, vstar[0]) - expr = Iu1 + f1 * y1 - assert extract_unique_domain(expr) == mesh1 \ No newline at end of file + u = Coefficient(V_1) + vstar = Argument(V_2.dual(), 0) + Iu = Interpolate(u, vstar) + assert extract_unique_domain(Iu) == mesh2 From c68448ad1021693c6d980dd8614d0072732f9d83 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 17 Jul 2025 14:22:12 +0100 Subject: [PATCH 23/35] fix imports --- ufl/action.py | 1 - ufl/integral.py | 1 + ufl/matrix.py | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/ufl/action.py b/ufl/action.py index 55dc611f2..4bd05cd66 100644 --- a/ufl/action.py +++ b/ufl/action.py @@ -9,7 +9,6 @@ from itertools import chain -from ufl import matrix # noqa 401 from ufl.algebra import Sum from ufl.argument import Argument, Coargument from ufl.coefficient import BaseCoefficient, Coefficient diff --git a/ufl/integral.py b/ufl/integral.py index 415060c96..1db3b72dc 100644 --- a/ufl/integral.py +++ b/ufl/integral.py @@ -32,6 +32,7 @@ class Integral: def __init__(self, integrand, integral_type, domain, subdomain_id, metadata, subdomain_data): """Initialise.""" + from ufl.domain import sort_domains if not isinstance(integrand, Expr): raise ValueError("Expecting integrand to be an Expr instance.") self._integrand = integrand diff --git a/ufl/matrix.py b/ufl/matrix.py index d548c0abd..990d5bcea 100644 --- a/ufl/matrix.py +++ b/ufl/matrix.py @@ -7,7 +7,6 @@ # # Modified by Nacime Bouziani, 2021-2022. -from ufl.argument import Argument from ufl.core.ufl_type import ufl_type from ufl.form import BaseForm from ufl.functionspace import AbstractFunctionSpace @@ -64,6 +63,7 @@ def ufl_function_spaces(self): def _analyze_form_arguments(self): """Define arguments of a matrix when considered as a form.""" + from ufl.argument import Argument self._arguments = ( Argument(self._ufl_function_spaces[0], 0), Argument(self._ufl_function_spaces[1], 1), From a55ebe759abbe36c96d75425f0d9d2f43e8ad638 Mon Sep 17 00:00:00 2001 From: Leo Date: Fri, 4 Jul 2025 12:28:02 +0100 Subject: [PATCH 24/35] cyclic import goose chase --- ufl/constant.py | 1 + 1 file changed, 1 insertion(+) diff --git a/ufl/constant.py b/ufl/constant.py index 384762550..539992ee4 100644 --- a/ufl/constant.py +++ b/ufl/constant.py @@ -20,6 +20,7 @@ class Constant(Terminal, Counted): def __init__(self, domain, shape=(), count=None): """Initalise.""" + from ufl.domain import as_domain Terminal.__init__(self) Counted.__init__(self, count, Constant) From 3c81eb74223311cccd78472c521a2f89b3507a54 Mon Sep 17 00:00:00 2001 From: Leo Date: Wed, 30 Jul 2025 13:45:21 +0100 Subject: [PATCH 25/35] add test for interpolate with mesh sequence --- test/test_unique_domain_extractor.py | 40 +++++++++++++++++++++++----- 1 file changed, 33 insertions(+), 7 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index 042f0df67..00c29673a 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -298,19 +298,45 @@ def test_extract_unique_domain_interpolate(): cell = triangle mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) - mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) - scalar_elem = LagrangeElement(cell, 1) vector_elem = LagrangeElement(cell, 1, (2,)) tensor_elem = LagrangeElement(cell, 1, (2, 2)) - V_1 = FunctionSpace(mesh1, scalar_elem) - V_2 = FunctionSpace(mesh2, vector_elem) - V_3 = FunctionSpace(mesh3, tensor_elem) + U = FunctionSpace(mesh1, vector_elem) + V = FunctionSpace(mesh2, tensor_elem) # # MixedFunctionSpace = V_3d x V_2d x V_1d # V = MixedFunctionSpace(V_1, V_2, V_3) - u = Coefficient(V_1) - vstar = Argument(V_2.dual(), 0) + u = Coefficient(U) + vstar = Argument(V.dual(), 0) Iu = Interpolate(u, vstar) assert extract_unique_domain(Iu) == mesh2 + + v = Coefficient(V) + v_test = Argument(V, 0) + u_test = Argument(U, 0) + form = inner(v_test, v) * Measure("ds", mesh2) + adjoint_I = Interpolate(u_test, form) # adjoint interpolation + assert extract_unique_domain(adjoint_I) == mesh1 + + +def test_extract_unique_domain_interpolate_meshsequence(): + cell = triangle + mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) + mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) + mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) + domain = MeshSequence([mesh1, mesh2, mesh3]) + + scalar_elem = LagrangeElement(cell, 1) + vector_elem = LagrangeElement(cell, 1, (2,)) + tensor_elem = LagrangeElement(cell, 1, (2, 2)) + mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) + + V = FunctionSpace(domain, mixed_elem) + + u1, u2, u3 = split(TestFunction(V)) + f1, f2, f3 = split(Coefficient(V)) + form = inner(u1, f1) * Measure("dx", mesh1) + + Iu = Interpolate(u2, form) # adjoint interpolation + assert extract_unique_domain(Iu) == mesh2 \ No newline at end of file From f6bd995633635fae5fc65aafb7dd4e1420595001 Mon Sep 17 00:00:00 2001 From: Leo Date: Wed, 30 Jul 2025 13:45:34 +0100 Subject: [PATCH 26/35] fix type check --- ufl/domain.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index c7d7aed14..3da94ab33 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -557,8 +557,9 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: # TODO: make this work for BaseForm # Action - Domain of 0th Argument in result # Leave AssembledMatrix and Adjoint for now - from ufl.form import Form, BaseForm + from ufl.form import Form from ufl.core.expr import Expr + from ufl.action import Action if isinstance(expr, Form): domains = set() @@ -573,7 +574,7 @@ def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: return domains.pop() else: raise ValueError(f"Form has multiple domains: {domains}") - elif isinstance(expr, BaseForm): + elif isinstance(expr, Action): return extract_unique_domain(expr.arguments()[0]) elif isinstance(expr, Expr): extractor = UniqueDomainExtractor() From c3933b3612f1f5e3e8d7175776ac9d5c4221fb50 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 28 Aug 2025 14:35:14 +0100 Subject: [PATCH 27/35] remove interpolate tests --- test/test_unique_domain_extractor.py | 48 ---------------------------- 1 file changed, 48 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index 00c29673a..eeb282177 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -292,51 +292,3 @@ def test_extract_unique_domain_action(): action_mixed = action(A_mixed, u_mixed) assert extract_unique_domain(action_mixed) == domain - - -def test_extract_unique_domain_interpolate(): - cell = triangle - mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) - mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) - vector_elem = LagrangeElement(cell, 1, (2,)) - tensor_elem = LagrangeElement(cell, 1, (2, 2)) - - U = FunctionSpace(mesh1, vector_elem) - V = FunctionSpace(mesh2, tensor_elem) - - # # MixedFunctionSpace = V_3d x V_2d x V_1d - # V = MixedFunctionSpace(V_1, V_2, V_3) - - u = Coefficient(U) - vstar = Argument(V.dual(), 0) - Iu = Interpolate(u, vstar) - assert extract_unique_domain(Iu) == mesh2 - - v = Coefficient(V) - v_test = Argument(V, 0) - u_test = Argument(U, 0) - form = inner(v_test, v) * Measure("ds", mesh2) - adjoint_I = Interpolate(u_test, form) # adjoint interpolation - assert extract_unique_domain(adjoint_I) == mesh1 - - -def test_extract_unique_domain_interpolate_meshsequence(): - cell = triangle - mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) - mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) - mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) - domain = MeshSequence([mesh1, mesh2, mesh3]) - - scalar_elem = LagrangeElement(cell, 1) - vector_elem = LagrangeElement(cell, 1, (2,)) - tensor_elem = LagrangeElement(cell, 1, (2, 2)) - mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) - - V = FunctionSpace(domain, mixed_elem) - - u1, u2, u3 = split(TestFunction(V)) - f1, f2, f3 = split(Coefficient(V)) - form = inner(u1, f1) * Measure("dx", mesh1) - - Iu = Interpolate(u2, form) # adjoint interpolation - assert extract_unique_domain(Iu) == mesh2 \ No newline at end of file From 283eb2fed6663470c009f09bdaf5193addbd06b5 Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 28 Aug 2025 14:48:31 +0100 Subject: [PATCH 28/35] tidy --- test/test_unique_domain_extractor.py | 5 ----- ufl/domain.py | 2 +- 2 files changed, 1 insertion(+), 6 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index eeb282177..c0516012d 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -8,13 +8,9 @@ MeshSequence, SpatialCoordinate, Measure, - Argument, - Interpolate, Constant, TrialFunction, - TrialFunctions, TestFunction, - TestFunctions, split, triangle, cos, @@ -23,7 +19,6 @@ FacetNormal, grad, action, - MixedFunctionSpace, Matrix ) diff --git a/ufl/domain.py b/ufl/domain.py index 3da94ab33..0ebf4f449 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -542,7 +542,7 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: return fs.ufl_domain() -def extract_unique_domain(expr: Union[Expr, Form]) -> AbstractDomain: +def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: """Extract the single unique domain from an expression. This works for expressions containing Indexed Arguments and Coefficients from From 71be2d4b50d1b6040d73b45c0669793cdf53731f Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 28 Aug 2025 15:01:29 +0100 Subject: [PATCH 29/35] remove interpolate case --- ufl/domain.py | 7 ------- 1 file changed, 7 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index 0ebf4f449..523464869 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -534,13 +534,6 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" ) - @process.register(Interpolate) - @DAGTraverser.postorder - def _(self, o: Expr, *operand_results) -> AbstractDomain: - """Process Interpolate.""" - fs = o.ufl_function_space() - return fs.ufl_domain() - def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: """Extract the single unique domain from an expression. From 4075d9ac79080c367f39f02cd1cae8a92dcc6aaf Mon Sep 17 00:00:00 2001 From: Leo Date: Thu, 28 Aug 2025 15:05:53 +0100 Subject: [PATCH 30/35] tidy --- ufl/domain.py | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index 523464869..36da73f39 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -433,7 +433,6 @@ def find_geometric_dimension(expr): class UniqueDomainExtractor(DAGTraverser): - # ExternalOperator or Interpolate: Domain is where we are mapping into. These are both BaseFormOperators. def __init__( self, compress: Union[bool, None] = True, @@ -464,10 +463,7 @@ def process(self, o: Expr) -> Expr: def _(self, o: Expr, *operand_results) -> AbstractDomain: """Process Indexed object by extracting the domain corresponding to the index.""" from ufl.functionspace import FunctionSpace - # Get the indexed expression and the multiindex expression, multiindex = o.ufl_operands - - # Get the domain from the first operand (the expression being indexed) expression_domain = operand_results[0] if isinstance(expression_domain, MeshSequence): @@ -475,23 +471,19 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: element = expression.ufl_element() if hasattr(element, 'sub_elements'): # Need to do this in case we have sub elements which are vector or tensor valued - offset = 0 + j = 0 for i, sub_element in enumerate(element.sub_elements): # Get the value size for this sub-element on its corresponding mesh sub_element_mesh = expression_domain.meshes[i] sub_element_fs = FunctionSpace(sub_element_mesh, sub_element) sub_element_size = sub_element_fs.value_size - if index < offset + sub_element_size: - # This index belongs to mesh i + if index < j + sub_element_size: return sub_element_mesh - offset += sub_element_size + j += sub_element_size raise ValueError(f"Index {index} out of range for mixed function space") else: - # If no sub elements we just grab the mesh return expression_domain.meshes[index] - - # If it's not a MeshSequence, just return the expression domain return expression_domain @process.register(Terminal) @@ -514,14 +506,11 @@ def _(self, o: Expr) -> AbstractDomain: @DAGTraverser.postorder def _(self, o: Expr, *operand_results) -> AbstractDomain: """Process Operator.""" - # Filter out None results (from operands that don't have domains) domains = [d for d in operand_results if d is not None] if not domains: - # No operands have domains return None elif len(domains) == 1: - # Only one operand has a domain return domains[0] else: # Multiple operands have domains - they should all be the same @@ -529,7 +518,6 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: if all(d == first_domain for d in domains): return first_domain else: - # If domains are not none and differ, then we don't have a unique domain raise ValueError( f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" ) @@ -550,6 +538,7 @@ def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: # TODO: make this work for BaseForm # Action - Domain of 0th Argument in result # Leave AssembledMatrix and Adjoint for now + # ExternalOperator or Interpolate: Domain is where we are mapping into. These are both BaseFormOperators. from ufl.form import Form from ufl.core.expr import Expr from ufl.action import Action From 05f620c21e0f85c0dcc9a0a75a202db11e8040e5 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 3 Sep 2025 15:21:13 +0100 Subject: [PATCH 31/35] `BaseForm` case --- test/test_unique_domain_extractor.py | 50 +++++++++++++++++----------- ufl/domain.py | 26 ++++++++++----- 2 files changed, 48 insertions(+), 28 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index c0516012d..a9b638d61 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -18,8 +18,10 @@ div, FacetNormal, grad, - action, - Matrix + Action, + Matrix, + Adjoint, + Interpolate ) from ufl.domain import extract_unique_domain @@ -259,31 +261,41 @@ def test_extract_unique_domain_repeated_meshes(): extract_unique_domain(u1 + u2) -def test_extract_unique_domain_action(): +def test_extract_unique_domain_baseform(): cell = triangle mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) mesh2 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=401) - mesh3 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=402) - domain = MeshSequence([mesh1, mesh2, mesh3]) scalar_elem = LagrangeElement(cell, 1) - vector_elem = LagrangeElement(cell, 1, (2,)) - tensor_elem = LagrangeElement(cell, 1, (2, 2)) - mixed_elem = MixedElement([scalar_elem, vector_elem, tensor_elem]) V1 = FunctionSpace(mesh1, scalar_elem) V2 = FunctionSpace(mesh2, scalar_elem) A = Matrix(V1, V2) - u = Coefficient(V2) - - action_Au = action(A, u) + assert extract_unique_domain(A) == mesh1 + + v = Coefficient(V2) + action_Au = Action(A, v) assert extract_unique_domain(action_Au) == mesh1 - - V_mixed1 = FunctionSpace(domain, mixed_elem) - V_mixed2 = FunctionSpace(domain, mixed_elem) - A_mixed = Matrix(V_mixed1, V_mixed2) - u_mixed = Coefficient(V_mixed2) - - action_mixed = action(A_mixed, u_mixed) - assert extract_unique_domain(action_mixed) == domain + Astar = Adjoint(A) + assert extract_unique_domain(Astar) == mesh2 + + v1 = TrialFunction(V1) + v2star = TestFunction(V2.dual()) + interp = Interpolate(v1, v2star) # V1 x V2^* -> R, equiv V1 -> V2 + assert extract_unique_domain(interp) == mesh2 + adjoint_interp = Adjoint(interp) # V2^* x V1 -> R, equiv V2^* -> V1^* + assert extract_unique_domain(adjoint_interp) == mesh1 + + cofunc = Coefficient(V2.dual()) + scalar = Action(cofunc, v) + assert extract_unique_domain(scalar) is None + + v = TestFunction(V2) + dx = Measure("dx", mesh2) + one_form = v * dx + formsum = cofunc + one_form + assert extract_unique_domain(formsum) is mesh2 + + two_form = interp * v * dx + assert extract_unique_domain(two_form) is mesh2 \ No newline at end of file diff --git a/ufl/domain.py b/ufl/domain.py index 36da73f39..3f9bf7e8e 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -10,7 +10,7 @@ import numbers from collections.abc import Iterable, Sequence -from typing import TYPE_CHECKING, Union +from typing import TYPE_CHECKING from functools import singledispatchmethod if TYPE_CHECKING: @@ -26,7 +26,7 @@ from ufl.core.operator import Operator from ufl.core.terminal import Terminal from ufl.indexed import Indexed -from ufl.core.interpolate import Interpolate +from ufl.core.base_form_operator import BaseFormOperator # Export list for ufl.classes __all_classes__ = ["AbstractDomain", "Mesh", "MeshView"] @@ -435,9 +435,9 @@ def find_geometric_dimension(expr): class UniqueDomainExtractor(DAGTraverser): def __init__( self, - compress: Union[bool, None] = True, - visited_cache: Union[dict[tuple, Expr], None] = None, - result_cache: Union[dict[Expr, Expr], None] = None, + compress: bool | None = True, + visited_cache: dict[tuple, Expr] | None = None, + result_cache: dict[Expr, Expr] | None = None, ) -> None: """Initialise.""" self._compress = compress @@ -487,6 +487,7 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: return expression_domain @process.register(Terminal) + @DAGTraverser.postorder def _(self, o: Expr) -> AbstractDomain: from ufl.coefficient import Coefficient from ufl.argument import Argument @@ -521,6 +522,11 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: raise ValueError( f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" ) + + @process.register(BaseFormOperator) + @DAGTraverser.postorder + def _(self, o: Expr, *operand_results) -> AbstractDomain: + pass def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: @@ -539,9 +545,8 @@ def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: # Action - Domain of 0th Argument in result # Leave AssembledMatrix and Adjoint for now # ExternalOperator or Interpolate: Domain is where we are mapping into. These are both BaseFormOperators. - from ufl.form import Form + from ufl.form import Form, BaseForm from ufl.core.expr import Expr - from ufl.action import Action if isinstance(expr, Form): domains = set() @@ -556,8 +561,11 @@ def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: return domains.pop() else: raise ValueError(f"Form has multiple domains: {domains}") - elif isinstance(expr, Action): - return extract_unique_domain(expr.arguments()[0]) + elif isinstance(expr, BaseForm): + if not expr.arguments(): + return None + else: + return extract_unique_domain(expr.arguments()[0]) elif isinstance(expr, Expr): extractor = UniqueDomainExtractor() return extractor(expr) From 43ff19090e03692bc5def5c28fad97ecd52b8ae5 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 3 Sep 2025 15:23:36 +0100 Subject: [PATCH 32/35] `BaseFormOperator` case --- ufl/domain.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/ufl/domain.py b/ufl/domain.py index 3f9bf7e8e..65dc4b7f1 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -526,7 +526,8 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: @process.register(BaseFormOperator) @DAGTraverser.postorder def _(self, o: Expr, *operand_results) -> AbstractDomain: - pass + fs = o.ufl_function_space() + return fs.ufl_domain() def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: @@ -541,10 +542,6 @@ def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: Returns: AbstractDomain: The unique domain extracted from the expression. """ - # TODO: make this work for BaseForm - # Action - Domain of 0th Argument in result - # Leave AssembledMatrix and Adjoint for now - # ExternalOperator or Interpolate: Domain is where we are mapping into. These are both BaseFormOperators. from ufl.form import Form, BaseForm from ufl.core.expr import Expr From 75bcbfad81cb6c6b64e2ce7f9a6336715fab6211 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Wed, 3 Sep 2025 15:26:49 +0100 Subject: [PATCH 33/35] lint --- test/test_unique_domain_extractor.py | 58 ++++++++++++++-------------- ufl/argument.py | 1 + ufl/checks.py | 1 + ufl/constant.py | 3 +- ufl/domain.py | 26 +++++++------ ufl/form.py | 5 ++- ufl/functionspace.py | 1 + ufl/integral.py | 1 - ufl/matrix.py | 1 + 9 files changed, 52 insertions(+), 45 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index a9b638d61..582332bd5 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -1,33 +1,33 @@ +import pytest from utils import FiniteElement, LagrangeElement, MixedElement -import pytest from ufl import ( + Action, + Adjoint, Coefficient, + Constant, + FacetNormal, FunctionSpace, + Interpolate, + Matrix, + Measure, Mesh, MeshSequence, SpatialCoordinate, - Measure, - Constant, - TrialFunction, TestFunction, - split, - triangle, + TrialFunction, cos, - inner, div, - FacetNormal, grad, - Action, - Matrix, - Adjoint, - Interpolate + inner, + split, + triangle, ) - from ufl.domain import extract_unique_domain from ufl.pullback import contravariant_piola, identity_pullback from ufl.sobolevspace import L2, HDiv + def test_extract_unique_domain(): cell = triangle elem0 = LagrangeElement(cell, 1) @@ -103,49 +103,49 @@ def test_extract_unique_domain_single_mesh(): """Test domain extraction for standard function spaces on a single mesh.""" cell = triangle mesh = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=200) - + # Test scalar elements P1 = LagrangeElement(cell, 1) V_scalar = FunctionSpace(mesh, P1) u_scalar = TrialFunction(V_scalar) f_scalar = Coefficient(V_scalar) - + assert extract_unique_domain(u_scalar) == mesh assert extract_unique_domain(f_scalar) == mesh - + P1_vec = LagrangeElement(cell, 1, (2,)) V_vector = FunctionSpace(mesh, P1_vec) u_vector = TrialFunction(V_vector) f_vector = Coefficient(V_vector) - + assert extract_unique_domain(u_vector) == mesh assert extract_unique_domain(f_vector) == mesh - + assert extract_unique_domain(u_vector[0]) == mesh assert extract_unique_domain(u_vector[1]) == mesh assert extract_unique_domain(f_vector[0]) == mesh assert extract_unique_domain(f_vector[1]) == mesh - + P1_tensor = LagrangeElement(cell, 1, (2, 2)) V_tensor = FunctionSpace(mesh, P1_tensor) u_tensor = TrialFunction(V_tensor) f_tensor = Coefficient(V_tensor) - + assert extract_unique_domain(u_tensor) == mesh assert extract_unique_domain(f_tensor) == mesh assert extract_unique_domain(u_tensor[0, 0]) == mesh assert extract_unique_domain(u_tensor[1, 1]) == mesh assert extract_unique_domain(f_tensor[0, 1]) == mesh - + x, y = SpatialCoordinate(mesh) expr1 = u_scalar + f_scalar expr2 = u_vector[0] + x expr3 = inner(u_vector, f_vector) - + assert extract_unique_domain(expr1) == mesh assert extract_unique_domain(expr2) == mesh assert extract_unique_domain(expr3) == mesh - + # Test forms dx = Measure("dx", mesh) form = inner(u_scalar, f_scalar) * dx @@ -153,7 +153,7 @@ def test_extract_unique_domain_single_mesh(): def test_extract_unique_domain_mixed_scalar_vector_tensor(): - """Test domain extraction for mixed function spaces + """Test domain extraction for mixed function spaces with scalar, vector, and tensor elements.""" cell = triangle mesh1 = Mesh(LagrangeElement(cell, 1, (2,)), ufl_id=400) @@ -177,16 +177,16 @@ def test_extract_unique_domain_mixed_scalar_vector_tensor(): assert extract_unique_domain(u_i) == domain[i] for i, f_i in enumerate((f_scalar, f_vector, f_tensor)): assert extract_unique_domain(f_i) == domain[i] - + for i in range(2): assert extract_unique_domain(u_vector[i]) == mesh2 assert extract_unique_domain(f_vector[i]) == mesh2 - + for i in range(2): for j in range(2): assert extract_unique_domain(u_tensor[i, j]) == mesh3 assert extract_unique_domain(f_tensor[i, j]) == mesh3 - + x1, y1 = SpatialCoordinate(mesh1) x2, y2 = SpatialCoordinate(mesh2) x3, y3 = SpatialCoordinate(mesh3) @@ -269,7 +269,7 @@ def test_extract_unique_domain_baseform(): V1 = FunctionSpace(mesh1, scalar_elem) V2 = FunctionSpace(mesh2, scalar_elem) - + A = Matrix(V1, V2) assert extract_unique_domain(A) == mesh1 @@ -298,4 +298,4 @@ def test_extract_unique_domain_baseform(): assert extract_unique_domain(formsum) is mesh2 two_form = interp * v * dx - assert extract_unique_domain(two_form) is mesh2 \ No newline at end of file + assert extract_unique_domain(two_form) is mesh2 diff --git a/ufl/argument.py b/ufl/argument.py index 01c9d1445..3f291bb47 100644 --- a/ufl/argument.py +++ b/ufl/argument.py @@ -285,6 +285,7 @@ def Arguments(function_space, number): Returns a tuple with the function components corresponding to the subelements. """ from ufl.split_functions import split + if isinstance(function_space, MixedFunctionSpace): return [ Argument(function_space.ufl_sub_space(i), number, i) diff --git a/ufl/checks.py b/ufl/checks.py index 9dab24615..d4c7a51df 100644 --- a/ufl/checks.py +++ b/ufl/checks.py @@ -40,6 +40,7 @@ def is_cellwise_constant(expr): def is_scalar_constant_expression(expr): """Check if an expression is a globally constant scalar expression.""" from ufl.geometry import GeometricQuantity + if is_python_scalar(expr): return True if expr.ufl_shape: diff --git a/ufl/constant.py b/ufl/constant.py index 539992ee4..c14e36a0a 100644 --- a/ufl/constant.py +++ b/ufl/constant.py @@ -8,8 +8,8 @@ from ufl.core.terminal import Terminal from ufl.core.ufl_type import ufl_type -from ufl.utils.counted import Counted from ufl.domain import as_domain +from ufl.utils.counted import Counted @ufl_type() @@ -21,6 +21,7 @@ class Constant(Terminal, Counted): def __init__(self, domain, shape=(), count=None): """Initalise.""" from ufl.domain import as_domain + Terminal.__init__(self) Counted.__init__(self, count, Constant) diff --git a/ufl/domain.py b/ufl/domain.py index 65dc4b7f1..f51012a00 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -10,23 +10,23 @@ import numbers from collections.abc import Iterable, Sequence -from typing import TYPE_CHECKING from functools import singledispatchmethod +from typing import TYPE_CHECKING if TYPE_CHECKING: from ufl.core.expr import Expr from ufl.finiteelement import AbstractFiniteElement # To avoid cyclic import when type-hinting. from ufl.form import Form from ufl.cell import AbstractCell +from ufl.core.base_form_operator import BaseFormOperator +from ufl.core.operator import Operator +from ufl.core.terminal import Terminal from ufl.core.ufl_id import attach_ufl_id from ufl.core.ufl_type import UFLObject -from ufl.corealg.traversal import traverse_unique_terminals -from ufl.sobolevspace import H1 from ufl.corealg.dag_traverser import DAGTraverser -from ufl.core.operator import Operator -from ufl.core.terminal import Terminal +from ufl.corealg.traversal import traverse_unique_terminals from ufl.indexed import Indexed -from ufl.core.base_form_operator import BaseFormOperator +from ufl.sobolevspace import H1 # Export list for ufl.classes __all_classes__ = ["AbstractDomain", "Mesh", "MeshView"] @@ -463,13 +463,14 @@ def process(self, o: Expr) -> Expr: def _(self, o: Expr, *operand_results) -> AbstractDomain: """Process Indexed object by extracting the domain corresponding to the index.""" from ufl.functionspace import FunctionSpace + expression, multiindex = o.ufl_operands expression_domain = operand_results[0] if isinstance(expression_domain, MeshSequence): index = multiindex[0]._value element = expression.ufl_element() - if hasattr(element, 'sub_elements'): + if hasattr(element, "sub_elements"): # Need to do this in case we have sub elements which are vector or tensor valued j = 0 for i, sub_element in enumerate(element.sub_elements): @@ -489,10 +490,11 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: @process.register(Terminal) @DAGTraverser.postorder def _(self, o: Expr) -> AbstractDomain: - from ufl.coefficient import Coefficient from ufl.argument import Argument - from ufl.geometry import GeometricQuantity + from ufl.coefficient import Coefficient from ufl.constant import Constant + from ufl.geometry import GeometricQuantity + if isinstance(o, (Coefficient, Argument)): fs = o.ufl_function_space() return fs.ufl_domain() @@ -522,12 +524,12 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: raise ValueError( f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" ) - + @process.register(BaseFormOperator) @DAGTraverser.postorder def _(self, o: Expr, *operand_results) -> AbstractDomain: fs = o.ufl_function_space() - return fs.ufl_domain() + return fs.ufl_domain() def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: @@ -542,8 +544,8 @@ def extract_unique_domain(expr: Expr | Form) -> AbstractDomain: Returns: AbstractDomain: The unique domain extracted from the expression. """ - from ufl.form import Form, BaseForm from ufl.core.expr import Expr + from ufl.form import BaseForm, Form if isinstance(expr, Form): domains = set() diff --git a/ufl/form.py b/ufl/form.py index 3f5ad0c17..ba2e07a74 100644 --- a/ufl/form.py +++ b/ufl/form.py @@ -17,7 +17,6 @@ from itertools import chain from ufl.checks import is_scalar_constant_expression - from ufl.constantvalue import Zero from ufl.core.expr import Expr, ufl_err_str from ufl.core.terminal import FormArgument @@ -40,6 +39,7 @@ def _sorted_integrals(integrals): stable signature computation. """ from ufl.domain import sort_domains + # Group integrals in multilevel dict by keys # [domain][integral_type][subdomain_id] integrals_dict = defaultdict(lambda: defaultdict(lambda: defaultdict(list))) @@ -400,6 +400,7 @@ def constants(self): def constant_numbering(self): """Return a contiguous numbering of constants in a mapping ``{constant:number}``.""" from ufl.constant import Constant + if self._constant_numbering is None: self._constant_numbering = { expr: num @@ -595,7 +596,7 @@ def __repr__(self): def _analyze_domains(self): """Analyze domains.""" - from ufl.domain import join_domains, sort_domains, extract_unique_domain + from ufl.domain import extract_unique_domain, join_domains, sort_domains # Collect integration domains. self._integration_domains = sort_domains( diff --git a/ufl/functionspace.py b/ufl/functionspace.py index 58616ee55..797d649b4 100644 --- a/ufl/functionspace.py +++ b/ufl/functionspace.py @@ -311,6 +311,7 @@ def ufl_element(self): def ufl_domains(self): """Return ufl domains.""" from ufl.domain import join_domains + domainlist = [] for s in self._ufl_function_spaces: domainlist.extend(s.ufl_domains()) diff --git a/ufl/integral.py b/ufl/integral.py index 1db3b72dc..415060c96 100644 --- a/ufl/integral.py +++ b/ufl/integral.py @@ -32,7 +32,6 @@ class Integral: def __init__(self, integrand, integral_type, domain, subdomain_id, metadata, subdomain_data): """Initialise.""" - from ufl.domain import sort_domains if not isinstance(integrand, Expr): raise ValueError("Expecting integrand to be an Expr instance.") self._integrand = integrand diff --git a/ufl/matrix.py b/ufl/matrix.py index 990d5bcea..42b009e58 100644 --- a/ufl/matrix.py +++ b/ufl/matrix.py @@ -64,6 +64,7 @@ def ufl_function_spaces(self): def _analyze_form_arguments(self): """Define arguments of a matrix when considered as a form.""" from ufl.argument import Argument + self._arguments = ( Argument(self._ufl_function_spaces[0], 0), Argument(self._ufl_function_spaces[1], 1), From e296ddc330dffdd8492d10f3f6a6bc8aede9513c Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 9 Sep 2025 11:32:43 +0100 Subject: [PATCH 34/35] tidy --- ufl/action.py | 1 + ufl/argument.py | 3 +-- ufl/constant.py | 2 -- ufl/corealg/dag_traverser.py | 8 +++++--- ufl/domain.py | 2 +- 5 files changed, 8 insertions(+), 8 deletions(-) diff --git a/ufl/action.py b/ufl/action.py index 4bd05cd66..55dc611f2 100644 --- a/ufl/action.py +++ b/ufl/action.py @@ -9,6 +9,7 @@ from itertools import chain +from ufl import matrix # noqa 401 from ufl.algebra import Sum from ufl.argument import Argument, Coargument from ufl.coefficient import BaseCoefficient, Coefficient diff --git a/ufl/argument.py b/ufl/argument.py index 3f291bb47..043f23dc8 100644 --- a/ufl/argument.py +++ b/ufl/argument.py @@ -21,6 +21,7 @@ from ufl.duals import is_dual, is_primal from ufl.form import BaseForm from ufl.functionspace import AbstractFunctionSpace, MixedFunctionSpace +from ufl.split_functions import split # Export list for ufl.classes (TODO: not actually classes: drop? these are in ufl.*) __all_classes__ = ["TestFunction", "TrialFunction", "TestFunctions", "TrialFunctions"] @@ -284,8 +285,6 @@ def Arguments(function_space, number): Returns a tuple with the function components corresponding to the subelements. """ - from ufl.split_functions import split - if isinstance(function_space, MixedFunctionSpace): return [ Argument(function_space.ufl_sub_space(i), number, i) diff --git a/ufl/constant.py b/ufl/constant.py index c14e36a0a..dfc65aa06 100644 --- a/ufl/constant.py +++ b/ufl/constant.py @@ -20,8 +20,6 @@ class Constant(Terminal, Counted): def __init__(self, domain, shape=(), count=None): """Initalise.""" - from ufl.domain import as_domain - Terminal.__init__(self) Counted.__init__(self, count, Constant) diff --git a/ufl/corealg/dag_traverser.py b/ufl/corealg/dag_traverser.py index 285e14e53..13ef83733 100644 --- a/ufl/corealg/dag_traverser.py +++ b/ufl/corealg/dag_traverser.py @@ -1,10 +1,12 @@ """Base class for dag traversers.""" +from __future__ import annotations from functools import singledispatchmethod, wraps -from typing import overload +from typing import overload, TYPE_CHECKING -from ufl.classes import Expr -from ufl.form import BaseForm +if TYPE_CHECKING: + from ufl.classes import Expr + from ufl.form import BaseForm class DAGTraverser: diff --git a/ufl/domain.py b/ufl/domain.py index f51012a00..758352274 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -522,7 +522,7 @@ def _(self, o: Expr, *operand_results) -> AbstractDomain: return first_domain else: raise ValueError( - f"Cannot extract unique domain from expression {o!r} with differing domains: {domains!r}" + f"Expression {o!r} has differing domains: {domains!r}" ) @process.register(BaseFormOperator) From 6a0ccd05b669d7345b577b8aa452e096833b5e69 Mon Sep 17 00:00:00 2001 From: Leo Collins Date: Tue, 9 Sep 2025 11:40:07 +0100 Subject: [PATCH 35/35] lint --- test/test_unique_domain_extractor.py | 2 +- ufl/corealg/dag_traverser.py | 3 ++- ufl/domain.py | 4 +++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/test/test_unique_domain_extractor.py b/test/test_unique_domain_extractor.py index 582332bd5..8d9c097c7 100644 --- a/test/test_unique_domain_extractor.py +++ b/test/test_unique_domain_extractor.py @@ -58,7 +58,7 @@ def test_extract_unique_domain(): assert extract_unique_domain(expr2) == mesh1 x2, y2 = SpatialCoordinate(mesh2) - with pytest.raises(ValueError) as e_info: + with pytest.raises(ValueError): _ = extract_unique_domain(u1 + u2) _ = extract_unique_domain(u1 + u2 + x2 * cos(x2 * u1)) diff --git a/ufl/corealg/dag_traverser.py b/ufl/corealg/dag_traverser.py index 13ef83733..f066afa4f 100644 --- a/ufl/corealg/dag_traverser.py +++ b/ufl/corealg/dag_traverser.py @@ -1,8 +1,9 @@ """Base class for dag traversers.""" from __future__ import annotations + from functools import singledispatchmethod, wraps -from typing import overload, TYPE_CHECKING +from typing import TYPE_CHECKING, overload if TYPE_CHECKING: from ufl.classes import Expr diff --git a/ufl/domain.py b/ufl/domain.py index 758352274..8514c4bbb 100644 --- a/ufl/domain.py +++ b/ufl/domain.py @@ -433,6 +433,8 @@ def find_geometric_dimension(expr): class UniqueDomainExtractor(DAGTraverser): + """Extract unique domain from an expression or BaseForm.""" + def __init__( self, compress: bool | None = True, @@ -495,7 +497,7 @@ def _(self, o: Expr) -> AbstractDomain: from ufl.constant import Constant from ufl.geometry import GeometricQuantity - if isinstance(o, (Coefficient, Argument)): + if isinstance(o, Coefficient | Argument): fs = o.ufl_function_space() return fs.ufl_domain() elif isinstance(o, GeometricQuantity):