From 0e3ad96b21bced84ab5ed7098cfaaab09c511c2a Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Fri, 31 Jul 2026 14:37:20 -0700
Subject: [PATCH 1/5] feat: generalise VectorField to N components

Support 3D vector fields via N-component VectorField, Component/Norm
operators, and matching discretisation/solution handling. Extracted from
the unstructured finite-volume work for review in isolation.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 CHANGELOG.md                                  |  4 +
 .../pybamm/discretisations/discretisation.py  | 72 +++++++++-----
 .../pybamm/expression_tree/unary_operators.py | 58 ++++++++++++
 .../pybamm/expression_tree/vector_field.py    | 94 ++++++++++++-------
 .../pybamm/src/pybamm/solvers/solution.py     | 34 +++++--
 packages/pybamm/tests/strategies/symbols.py   | 15 +++
 .../test_tensor_field.py                      | 77 +++++++++++++++
 7 files changed, 287 insertions(+), 67 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 79dcab15e8..212d4ea8e0 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Features
+
+- Generalised `VectorField` to N components and added `Component`/`Norm` operators for multi-dimensional vector fields.
+
 ## Bug fixes
 
 - `BatchStudy.solve` no longer ignores its `solver` argument: previously the loop over study inputs shadowed it, so a caller-supplied solver was silently dropped. A solver from `BatchStudy(solvers=...)` still takes precedence. ([#5677](https://github.com/pybamm-team/PyBaMM/pull/5677))
diff --git a/packages/pybamm/src/pybamm/discretisations/discretisation.py b/packages/pybamm/src/pybamm/discretisations/discretisation.py
index 779b64c2a2..86cacc0dd4 100644
--- a/packages/pybamm/src/pybamm/discretisations/discretisation.py
+++ b/packages/pybamm/src/pybamm/discretisations/discretisation.py
@@ -882,7 +882,11 @@ def process_symbol(self, symbol):
 
         # Assign mesh as an attribute to the processed variable
         if symbol.domain != []:
-            discretised_symbol.mesh = self.mesh[symbol.domain]
+            mesh_for_symbol = self.mesh[symbol.domain]
+            discretised_symbol.mesh = mesh_for_symbol
+            if isinstance(discretised_symbol, pybamm.VectorField):
+                for comp in discretised_symbol._components:
+                    comp.mesh = mesh_for_symbol
         else:
             discretised_symbol.mesh = None
 
@@ -934,23 +938,31 @@ def _process_symbol(self, symbol):
                 if isinstance(disc_left, pybamm.VectorField) or isinstance(
                     disc_right, pybamm.VectorField
                 ):
+                    if isinstance(disc_left, pybamm.VectorField):
+                        n = disc_left.n_components
+                    else:
+                        n = disc_right.n_components
                     if not isinstance(disc_right, pybamm.VectorField):
-                        disc_right = pybamm.VectorField(disc_right, disc_right)
+                        disc_right = pybamm.VectorField(*[disc_right] * n)
                     if not isinstance(disc_left, pybamm.VectorField):
-                        disc_left = pybamm.VectorField(disc_left, disc_left)
-                    else:  # both are vector fields already
-                        pass
-                    disc_lr = pybamm.simplify_if_constant(
-                        symbol.create_copy(
-                            new_children=[disc_left.lr_field, disc_right.lr_field]
-                        )
-                    )
-                    disc_tb = pybamm.simplify_if_constant(
-                        symbol.create_copy(
-                            new_children=[disc_left.tb_field, disc_right.tb_field]
+                        disc_left = pybamm.VectorField(*[disc_left] * n)
+                    new_comps = [
+                        pybamm.simplify_if_constant(
+                            symbol.create_copy(
+                                new_children=[
+                                    disc_left._components[k],
+                                    disc_right._components[k],
+                                ]
+                            )
                         )
-                    )
-                    return pybamm.VectorField(disc_lr, disc_tb)
+                        for k in range(n)
+                    ]
+                    result = pybamm.VectorField(*new_comps)
+                    for src in (disc_left, disc_right):
+                        if hasattr(src, "_disc_state_vector"):
+                            result._disc_state_vector = src._disc_state_vector
+                            break
+                    return result
 
                 return pybamm.simplify_if_constant(
                     symbol.create_copy(new_children=[disc_left, disc_right])
@@ -1092,6 +1104,18 @@ def _process_symbol(self, symbol):
             elif isinstance(symbol, pybamm.NotConstant):
                 # After discretisation, we can make the symbol constant
                 return disc_child
+            elif isinstance(symbol, pybamm.Component):
+                if not isinstance(disc_child, pybamm.VectorField):
+                    raise ValueError("Component can only be applied to a VectorField")
+                return disc_child._components[symbol.index]
+            elif isinstance(symbol, pybamm.Norm):
+                if not isinstance(disc_child, pybamm.VectorField):
+                    raise ValueError("Norm can only be applied to a VectorField")
+                result = None
+                for comp in disc_child._components:
+                    sq = comp**2
+                    result = sq if result is None else result + sq
+                return result**0.5
             elif isinstance(symbol, pybamm.Magnitude):
                 if not isinstance(disc_child, pybamm.VectorField):
                     raise ValueError("Magnitude can only be applied to a vector field")
@@ -1104,10 +1128,14 @@ def _process_symbol(self, symbol):
                     raise ValueError("Invalid direction")
             else:
                 if isinstance(disc_child, pybamm.VectorField):
-                    return pybamm.VectorField(
-                        symbol.create_copy(new_children=[disc_child.lr_field]),
-                        symbol.create_copy(new_children=[disc_child.tb_field]),
-                    )
+                    new_comps = [
+                        symbol.create_copy(new_children=[c])
+                        for c in disc_child._components
+                    ]
+                    result = pybamm.VectorField(*new_comps)
+                    if hasattr(disc_child, "_disc_state_vector"):
+                        result._disc_state_vector = disc_child._disc_state_vector
+                    return result
                 else:
                     return symbol.create_copy(new_children=[disc_child])
 
@@ -1181,10 +1209,8 @@ def _process_symbol(self, symbol):
             )
 
         elif isinstance(symbol, pybamm.VectorField):
-            # VectorField is a subclass of TensorField, handle it first for specificity
-            left_symbol = self.process_symbol(symbol.lr_field)
-            right_symbol = self.process_symbol(symbol.tb_field)
-            return symbol.create_copy(new_children=[left_symbol, right_symbol])
+            processed = [self.process_symbol(c) for c in symbol._components]
+            return symbol.create_copy(new_children=processed)
 
         elif isinstance(symbol, pybamm.TensorField):
             # General TensorField handling (rank-2 tensors)
diff --git a/packages/pybamm/src/pybamm/expression_tree/unary_operators.py b/packages/pybamm/src/pybamm/expression_tree/unary_operators.py
index 5e509d0684..be3f2df585 100644
--- a/packages/pybamm/src/pybamm/expression_tree/unary_operators.py
+++ b/packages/pybamm/src/pybamm/expression_tree/unary_operators.py
@@ -1509,6 +1509,54 @@ def _unary_new_copy(self, child, perform_simplifications=True):
         return self.__class__(child, self.direction)
 
 
+class Component(UnaryOperator):
+    """
+    Extract component *index* from a VectorField.
+
+    Parameters
+    ----------
+    child : :class:`pybamm.Symbol`
+        A VectorField symbol.
+    index : int
+        Zero-based component index.
+    """
+
+    def __init__(self, child, index):
+        super().__init__(f"component({index})", child)
+        self.index = index
+
+    def to_json(self):
+        return {
+            "name": self.name,
+            "domains": self.domains,
+            "index": self.index,
+        }
+
+    @classmethod
+    def _from_json(cls, snippet):
+        return cls(snippet["children"][0], snippet["index"])
+
+    def _unary_new_copy(self, child, perform_simplifications=True):
+        return self.__class__(child, self.index)
+
+
+class Norm(UnaryOperator):
+    """
+    Euclidean norm of a VectorField: ``sqrt(sum(comp_i ** 2))``.
+
+    Parameters
+    ----------
+    child : :class:`pybamm.Symbol`
+        A VectorField symbol.
+    """
+
+    def __init__(self, child):
+        super().__init__("norm", child)
+
+    def _unary_new_copy(self, child, perform_simplifications=True):
+        return self.__class__(child)
+
+
 class Upwind(UpwindDownwind):
     """
     Upwinding operator. To be used if flow velocity is positive (left to right).
@@ -1763,6 +1811,16 @@ def sign(symbol):
     return pybamm.simplify_if_constant(Sign(symbol))
 
 
+def component(symbol, index):
+    """Convenience function for creating a :class:`Component`."""
+    return Component(symbol, index)
+
+
+def norm(symbol):
+    """Convenience function for creating a :class:`Norm`."""
+    return Norm(symbol)
+
+
 def smooth_absolute_value(symbol, k):
     """
     Smooth approximation to the absolute value function. k is the smoothing parameter,
diff --git a/packages/pybamm/src/pybamm/expression_tree/vector_field.py b/packages/pybamm/src/pybamm/expression_tree/vector_field.py
index 031674d99e..4004c4d47c 100644
--- a/packages/pybamm/src/pybamm/expression_tree/vector_field.py
+++ b/packages/pybamm/src/pybamm/expression_tree/vector_field.py
@@ -1,9 +1,11 @@
 """
-VectorField class - a rank-1 tensor field for 2D simulations.
+VectorField class - a rank-1 tensor field with N components.
 """
 
 from __future__ import annotations
 
+import casadi
+
 import pybamm
 from pybamm.expression_tree.tensor_field import TensorField
 
@@ -13,22 +15,29 @@ class VectorField(TensorField):
     A node in the expression tree representing a vector field.
 
     VectorField is a convenience subclass of TensorField for rank-1 tensors
-    with two components (lr and tb directions in 2D).
+    with N >= 2 components.  Components are stored by integer index; the
+    properties ``lr_field``, ``tb_field``, and ``fb_field`` are backward-
+    compatible aliases for ``[0]``, ``[1]``, and ``[2]``.
 
     Parameters
     ----------
-    lr_field : pybamm.Symbol
-        The left-right (x) component of the vector field.
-    tb_field : pybamm.Symbol
-        The top-bottom (y) component of the vector field.
+    *components : pybamm.Symbol
+        Two or more component symbols, all sharing the same domain.
     """
 
-    def __init__(self, lr_field, tb_field):
-        if lr_field.domain != tb_field.domain:
-            raise ValueError("lr_field and tb_field must have the same domain")
-        # Initialize as a rank-1 TensorField with two components
-        super().__init__([lr_field, tb_field], domain=lr_field.domain)
-        # Override the name to maintain backward compatibility
+    def __init__(self, *components):
+        if len(components) < 2:
+            raise ValueError(
+                f"VectorField requires at least 2 components, got {len(components)}"
+            )
+        ref_domain = components[0].domain
+        for i, c in enumerate(components[1:], start=1):
+            if c.domain != ref_domain:
+                raise ValueError(
+                    f"All components must have the same domain: "
+                    f"component {i} has {c.domain}, expected {ref_domain}"
+                )
+        super().__init__(list(components), domain=ref_domain)
         self.name = "vector_field"
 
     @classmethod
@@ -36,44 +45,61 @@ def _from_json(cls, snippet):
         # Two positional args, not a single list -- override TensorField._from_json.
         return cls(snippet["children"][0], snippet["children"][1])
 
+    @property
+    def n_components(self):
+        """Number of vector components."""
+        return len(self._components)
+
+    # ---- backward-compatible aliases for structured-grid directions ----
+
     @property
     def lr_field(self):
-        """The left-right (x) component of the vector field."""
+        """Component 0 (left-right / x)."""
         return self._components[0]
 
     @property
     def tb_field(self):
-        """The top-bottom (y) component of the vector field."""
+        """Component 1 (top-bottom / y)."""
         return self._components[1]
 
+    @property
+    def fb_field(self):
+        """Component 2 (front-back / z).  Only valid for 3-component fields."""
+        if len(self._components) < 3:
+            raise AttributeError(
+                "fb_field requires at least 3 components; this VectorField has "
+                f"{len(self._components)}"
+            )
+        return self._components[2]
+
     def create_copy(
         self,
         new_children: list[pybamm.Symbol] | None = None,
         perform_simplifications: bool = True,
     ):
-        """Create a copy of this vector field with optional new children."""
         if new_children is None:
             new_children = [
-                self.lr_field.create_copy(
-                    perform_simplifications=perform_simplifications
-                ),
-                self.tb_field.create_copy(
-                    perform_simplifications=perform_simplifications
-                ),
+                c.create_copy(perform_simplifications=perform_simplifications)
+                for c in self._components
             ]
         return VectorField(*new_children)
 
+    def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
+        """See :meth:`pybamm.Symbol._to_casadi()`."""
+        return casadi.vertcat(
+            *[
+                c._to_casadi_inner(t, y, y_dot, inputs, casadi_symbols)
+                for c in self._components
+            ]
+        )
+
     def evaluates_on_edges(self, dimension: str) -> bool:
-        """Check if components evaluate on edges.
-
-        Overrides TensorField to provide more specific error message.
-        """
-        left_evaluates_on_edges = self.lr_field.evaluates_on_edges(dimension)
-        right_evaluates_on_edges = self.tb_field.evaluates_on_edges(dimension)
-        if left_evaluates_on_edges == right_evaluates_on_edges:
-            return left_evaluates_on_edges
-        else:
-            raise ValueError(
-                "lr_field and tb_field must either both evaluate on edges "
-                "or both not evaluate on edges"
-            )
+        statuses = [c.evaluates_on_edges(dimension) for c in self._components]
+        if all(statuses):
+            return True
+        if not any(statuses):
+            return False
+        raise ValueError(
+            "All VectorField components must either all evaluate on edges "
+            "or none evaluate on edges"
+        )
diff --git a/packages/pybamm/src/pybamm/solvers/solution.py b/packages/pybamm/src/pybamm/solvers/solution.py
index a89af76cb2..4e43b0390b 100644
--- a/packages/pybamm/src/pybamm/solvers/solution.py
+++ b/packages/pybamm/src/pybamm/solvers/solution.py
@@ -748,16 +748,30 @@ def _update_variable(self, name: str):
                     "solve. Please re-run the solve with `output_variables` set to "
                     "include this variable."
                 )
-            var_casadi, var_pybamm, time_integral = self._update_model_variable(
-                model,
-                _var_pybamm,
-                inputs=inputs,
-                ys_shape=ys.shape,
-                time_integral=time_integral,
-                cache_key=name,
-            )
-            vars_pybamm[i] = var_pybamm
-            vars_casadi[i] = var_casadi
+            if isinstance(_var_pybamm, pybamm.VectorField):
+                comp_casadi = []
+                for k, comp in enumerate(_var_pybamm._components):
+                    cc, _, _ = self._update_model_variable(
+                        model,
+                        comp,
+                        inputs=inputs,
+                        ys_shape=ys.shape,
+                        time_integral=None,
+                        cache_key=f"{name}[{k}]",
+                    )
+                    comp_casadi.append(cc)
+                vars_casadi[i] = comp_casadi
+            else:
+                var_casadi, var_pybamm, time_integral = self._update_model_variable(
+                    model,
+                    _var_pybamm,
+                    inputs=inputs,
+                    ys_shape=ys.shape,
+                    time_integral=time_integral,
+                    cache_key=name,
+                )
+                vars_pybamm[i] = var_pybamm
+                vars_casadi[i] = var_casadi
         var = pybamm.process_variable(
             name, vars_pybamm, vars_casadi, self, time_integral=time_integral
         )
diff --git a/packages/pybamm/tests/strategies/symbols.py b/packages/pybamm/tests/strategies/symbols.py
index 270519e7c1..33ab2751ab 100644
--- a/packages/pybamm/tests/strategies/symbols.py
+++ b/packages/pybamm/tests/strategies/symbols.py
@@ -670,6 +670,18 @@ def _magnitude_branch(
     )
 
 
+
+def _component_branch(
+    _child_strategy: st.SearchStrategy[pybamm.Symbol],
+) -> st.SearchStrategy[pybamm.Component]:
+    """Component(child, index) — domain-bearing child, zero-based component index."""
+    return st.builds(
+        pybamm.Component,
+        _any_domain_leaves(),
+        st.integers(min_value=0, max_value=2),
+    )
+
+
 def _discrete_time_sum_branch(
     _child_strategy: st.SearchStrategy[pybamm.Symbol],
 ) -> st.SearchStrategy[pybamm.DiscreteTimeSum]:
@@ -981,6 +993,9 @@ def _vector_branch(
         pybamm.UpwindDownwind2D: _upwind_downwind_2d_branch,
         pybamm.NodeToEdge2D: _node_to_edge_2d_branch,
         pybamm.Magnitude: _magnitude_branch,
+        pybamm.Component: _component_branch,
+        # Norm: (self, child) only — round-trips via the generic unary hook.
+        pybamm.Norm: lambda _children: _any_domain_leaves().map(pybamm.Norm),
         pybamm.DiscreteTimeData: _discrete_time_data_branch,
         pybamm.DiscreteTimeSum: _discrete_time_sum_branch,
         pybamm.SizeAverage: _size_average_branch,
diff --git a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
index 1568b20561..8201835925 100644
--- a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
+++ b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
@@ -126,6 +126,46 @@ def test_rank2_evaluates_on_edges(self):
         t = TensorField([[a, b], [c, d]])
         assert t.evaluates_on_edges("primary") is False
 
+    def test_components_property(self):
+        """Accessing the components property returns nested structure."""
+        a, b = pybamm.Scalar(1), pybamm.Scalar(2)
+        t = TensorField([a, b])
+        assert t.components == [a, b]
+
+    def test_rank1_tuple_index(self):
+        """Rank-1 tensor accepts single-element tuple index."""
+        a, b = pybamm.Scalar(1), pybamm.Scalar(2)
+        t = TensorField([a, b])
+        assert t[(0,)] == a
+
+    def test_rank1_too_many_indices_raises(self):
+        """Rank-1 tensor raises for multi-element tuple index."""
+        a, b = pybamm.Scalar(1), pybamm.Scalar(2)
+        t = TensorField([a, b])
+        with pytest.raises(IndexError, match="Too many indices for rank-1"):
+            t[(0, 1)]
+
+    def test_rank2_single_element_tuple_returns_row(self):
+        """Rank-2 tensor with single-element tuple returns row."""
+        a, b, c, d = [pybamm.Scalar(i) for i in range(4)]
+        t = TensorField([[a, b], [c, d]])
+        assert t[(0,)] == [a, b]
+
+    def test_rank2_too_many_indices_raises(self):
+        """Rank-2 tensor raises for 3+ element tuple index."""
+        a, b, c, d = [pybamm.Scalar(i) for i in range(4)]
+        t = TensorField([[a, b], [c, d]])
+        with pytest.raises(IndexError, match="Too many indices for rank-2"):
+            t[(0, 1, 2)]
+
+    def test_rank2_evaluates_on_edges_all_true(self):
+        """Rank-2 evaluates_on_edges returns True when all components are on edges."""
+        a, b, c, d = [pybamm.Scalar(i) for i in range(4)]
+        t = TensorField([[a, b], [c, d]])
+        for child in t.children:
+            child._evaluates_on_edges = lambda _: True
+        assert t.evaluates_on_edges("primary") is True
+
 
 class TestVectorFieldInheritance:
     """Tests for VectorField inheriting from TensorField."""
@@ -156,6 +196,43 @@ def test_vectorfield_domain_validation(self):
         with pytest.raises(ValueError, match="same domain"):
             pybamm.VectorField(a, b)
 
+    def test_vectorfield_requires_two_components(self):
+        """VectorField with fewer than 2 components raises."""
+        with pytest.raises(ValueError, match="requires at least 2 components"):
+            pybamm.VectorField(pybamm.Scalar(1))
+
+    def test_vectorfield_fb_field_three_components(self):
+        """fb_field returns 3rd component for 3-component VectorField."""
+        a, b, c = pybamm.Scalar(1), pybamm.Scalar(2), pybamm.Scalar(3)
+        vf = pybamm.VectorField(a, b, c)
+        assert vf.fb_field == c
+        assert vf.n_components == 3
+
+    def test_vectorfield_fb_field_raises_when_missing(self):
+        """fb_field on 2-component VectorField raises AttributeError."""
+        vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
+        with pytest.raises(AttributeError, match="fb_field requires at least 3"):
+            _ = vf.fb_field
+
+    def test_vectorfield_evaluates_on_edges_all_true(self):
+        """VectorField evaluates_on_edges returns True when all on edges."""
+        vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
+        vf.lr_field._evaluates_on_edges = lambda _: True
+        vf.tb_field._evaluates_on_edges = lambda _: True
+        assert vf.evaluates_on_edges("primary") is True
+
+    def test_vectorfield_to_casadi(self):
+        """VectorField _to_casadi concatenates components via vertcat."""
+        import casadi
+
+        a, b = pybamm.Scalar(1.0), pybamm.Scalar(2.0)
+        vf = pybamm.VectorField(a, b)
+        mx = vf.to_casadi()
+        assert isinstance(mx, casadi.MX)
+        f = casadi.Function("f", [], [mx])
+        out = f.call([])
+        np.testing.assert_array_equal(np.array(out[0]).flatten(), [1.0, 2.0])
+
 
 class TestTensorProduct:
     """Tests for TensorProduct operator."""

From 1414f81ca5462c919b3ed65def45e50a1a189081 Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Fri, 31 Jul 2026 21:39:53 +0000
Subject: [PATCH 2/5] style: pre-commit fixes

---
 packages/pybamm/tests/strategies/symbols.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/packages/pybamm/tests/strategies/symbols.py b/packages/pybamm/tests/strategies/symbols.py
index 33ab2751ab..4ad7a89627 100644
--- a/packages/pybamm/tests/strategies/symbols.py
+++ b/packages/pybamm/tests/strategies/symbols.py
@@ -670,7 +670,6 @@ def _magnitude_branch(
     )
 
 
-
 def _component_branch(
     _child_strategy: st.SearchStrategy[pybamm.Symbol],
 ) -> st.SearchStrategy[pybamm.Component]:

From 0e1522c6a763b1f4b235ef27baca6a2cb6bf5f86 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Fri, 31 Jul 2026 16:38:50 -0700
Subject: [PATCH 3/5] test: cover Component/Norm and VectorField solution paths

Add unit tests for the uncovered patch lines: Component/Norm
discretisation (success and error paths), _unary_new_copy and the
component/norm convenience functions, _disc_state_vector propagation
through binary/unary operators, and the per-component casadi handling
for VectorField variables in Solution.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
---
 .../test_vector_field_and_magnitude.py        | 36 +++++++++++++++++++
 .../tests/unit/test_solvers/test_solution.py  | 22 ++++++++++++
 .../test_tensor_field.py                      | 24 +++++++++++++
 3 files changed, 82 insertions(+)

diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py b/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py
index f1f798775a..696f36abf0 100644
--- a/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py
+++ b/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py
@@ -64,3 +64,39 @@ def test_vector_field_and_magnitude(self, mesh_2d):
 
         with pytest.raises(ValueError, match=r"Invalid direction"):
             disc.process_symbol(pybamm.Magnitude(vector_field, "asdf"))
+
+    def test_component_and_norm_discretisation(self, mesh_2d):
+        spatial_methods = {"macroscale": pybamm.FiniteVolume2D()}
+        disc = pybamm.Discretisation(mesh_2d, spatial_methods)
+        vector_field = pybamm.VectorField(pybamm.Scalar(3), pybamm.Scalar(4))
+
+        comp_0 = disc.process_symbol(pybamm.component(vector_field, 0))
+        comp_1 = disc.process_symbol(pybamm.component(vector_field, 1))
+        assert comp_0.evaluate() == 3
+        assert comp_1.evaluate() == 4
+
+        norm = disc.process_symbol(pybamm.norm(vector_field))
+        assert norm.evaluate() == pytest.approx(5.0)
+
+        with pytest.raises(ValueError, match=r"Component can only be applied"):
+            disc.process_symbol(pybamm.component(pybamm.Scalar(1), 0))
+
+        with pytest.raises(ValueError, match=r"Norm can only be applied"):
+            disc.process_symbol(pybamm.norm(pybamm.Scalar(1)))
+
+    def test_disc_state_vector_propagation(self, mesh_2d):
+        # binary and unary operators on a discretised VectorField must carry
+        # the _disc_state_vector attribute over to the result
+        spatial_methods = {"macroscale": pybamm.FiniteVolume2D()}
+        disc = pybamm.Discretisation(mesh_2d, spatial_methods)
+        vector_field = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
+        disc_vf = disc.process_symbol(vector_field)
+        marker = pybamm.StateVector(slice(0, 1))
+        disc_vf._disc_state_vector = marker
+
+        one = pybamm.Constant(1, "one")
+        disc_sum = disc.process_symbol(vector_field + one)
+        assert disc_sum._disc_state_vector is marker
+
+        disc_neg = disc.process_symbol(-vector_field)
+        assert disc_neg._disc_state_vector is marker
diff --git a/packages/pybamm/tests/unit/test_solvers/test_solution.py b/packages/pybamm/tests/unit/test_solvers/test_solution.py
index 95cf1597d8..367485c265 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_solution.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_solution.py
@@ -512,6 +512,28 @@ def test_observable_computed_lazily(self):
         assert sol_sum.observable == expected  # computed on access
         assert sol_sum._observable == expected  # and cached
 
+    def test_vector_field_variable(self, mesh_2d):
+        # a VectorField model variable is processed component-by-component,
+        # storing a list of casadi functions for the downstream processed
+        # variable classes to consume
+        model = pybamm.BaseModel()
+        var = pybamm.Variable(
+            "var", domain=["negative electrode", "separator", "positive electrode"]
+        )
+        model.rhs = {var: pybamm.Scalar(-1)}
+        model.initial_conditions = {var: pybamm.Scalar(1)}
+        model.variables = {"var": var, "flux": pybamm.VectorField(var, 2 * var)}
+
+        disc = pybamm.Discretisation(mesh_2d, {"macroscale": pybamm.FiniteVolume2D()})
+        disc.process_model(model)
+        solution = pybamm.IDAKLUSolver().solve(model, np.linspace(0, 1, 5))
+
+        flux = solution["flux"]
+        assert isinstance(flux.base_variables[0], pybamm.VectorField)
+        casadi_components = flux.base_variables_casadi[0]
+        assert isinstance(casadi_components, list)
+        assert len(casadi_components) == 2
+
     def test_add_solutions_different_models(self):
         # Set up first solution
         t1 = np.linspace(0, 1)
diff --git a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
index 8201835925..604ea629a5 100644
--- a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
+++ b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
@@ -234,6 +234,30 @@ def test_vectorfield_to_casadi(self):
         np.testing.assert_array_equal(np.array(out[0]).flatten(), [1.0, 2.0])
 
 
+class TestComponentAndNorm:
+    """Tests for the Component and Norm operators on VectorFields."""
+
+    def test_component_convenience_function_and_copy(self):
+        """pybamm.component creates a Component; copying preserves the index."""
+        vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
+        comp = pybamm.component(vf, 1)
+        assert isinstance(comp, pybamm.Component)
+        assert comp.index == 1
+        copy = comp.create_copy()
+        assert isinstance(copy, pybamm.Component)
+        assert copy.index == 1
+        assert copy == comp
+
+    def test_norm_convenience_function_and_copy(self):
+        """pybamm.norm creates a Norm; copying preserves structure."""
+        vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
+        norm = pybamm.norm(vf)
+        assert isinstance(norm, pybamm.Norm)
+        copy = norm.create_copy()
+        assert isinstance(copy, pybamm.Norm)
+        assert copy == norm
+
+
 class TestTensorProduct:
     """Tests for TensorProduct operator."""
 

From e347d972931dafc8930e2b3b86794e9d21e87392 Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Mon, 3 Aug 2026 11:54:00 -0700
Subject: [PATCH 4/5] Address review feedback on VectorField N-comp
 generalisation.

Guard mismatched component counts, raise a clear error for structured-2D
VectorField solution reads, and drop unused fb_field / convenience helpers.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 CHANGELOG.md                                  |  2 +-
 .../api/expression_tree/unary_operator.rst    |  6 ++
 .../pybamm/discretisations/discretisation.py  | 84 ++++++++++---------
 .../pybamm/expression_tree/unary_operators.py | 10 ---
 .../pybamm/expression_tree/vector_field.py    | 33 +++-----
 .../src/pybamm/solvers/processed_variable.py  |  6 ++
 .../pybamm/src/pybamm/solvers/solution.py     |  2 +-
 .../test_vector_field_and_magnitude.py        | 28 +++++--
 .../tests/unit/test_solvers/test_solution.py  | 13 ++-
 .../test_tensor_field.py                      | 24 ++----
 10 files changed, 106 insertions(+), 102 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 212d4ea8e0..88407fd6ff 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,7 +2,7 @@
 
 ## Features
 
-- Generalised `VectorField` to N components and added `Component`/`Norm` operators for multi-dimensional vector fields.
+- Generalised `VectorField` to N components and added `Component`/`Norm` operators for multi-dimensional vector fields. ([#5686](https://github.com/pybamm-team/PyBaMM/pull/5686))
 
 ## Bug fixes
 
diff --git a/docs/source/api/expression_tree/unary_operator.rst b/docs/source/api/expression_tree/unary_operator.rst
index e6a3cbe554..38dc345bdd 100644
--- a/docs/source/api/expression_tree/unary_operator.rst
+++ b/docs/source/api/expression_tree/unary_operator.rst
@@ -73,6 +73,12 @@ Unary Operators
 .. autoclass:: pybamm.Downwind
   :members:
 
+.. autoclass:: pybamm.Component
+  :members:
+
+.. autoclass:: pybamm.Norm
+  :members:
+
 .. autofunction:: pybamm.grad
 
 .. autofunction:: pybamm.div
diff --git a/packages/pybamm/src/pybamm/discretisations/discretisation.py b/packages/pybamm/src/pybamm/discretisations/discretisation.py
index 86cacc0dd4..046170af6e 100644
--- a/packages/pybamm/src/pybamm/discretisations/discretisation.py
+++ b/packages/pybamm/src/pybamm/discretisations/discretisation.py
@@ -885,7 +885,7 @@ def process_symbol(self, symbol):
             mesh_for_symbol = self.mesh[symbol.domain]
             discretised_symbol.mesh = mesh_for_symbol
             if isinstance(discretised_symbol, pybamm.VectorField):
-                for comp in discretised_symbol._components:
+                for comp in discretised_symbol.components:
                     comp.mesh = mesh_for_symbol
         else:
             discretised_symbol.mesh = None
@@ -904,6 +904,37 @@ def process_symbol(self, symbol):
 
         return discretised_symbol
 
+    def _process_vector_field_binary(self, symbol, disc_left, disc_right):
+        """Broadcast a scalar side, then apply ``symbol`` component-wise."""
+        left_is_vf = isinstance(disc_left, pybamm.VectorField)
+        right_is_vf = isinstance(disc_right, pybamm.VectorField)
+        if left_is_vf and right_is_vf:
+            if disc_left.n_components != disc_right.n_components:
+                raise pybamm.DiscretisationError(
+                    f"Cannot combine VectorFields with {disc_left.n_components} and "
+                    f"{disc_right.n_components} components"
+                )
+            n = disc_left.n_components
+        elif left_is_vf:
+            n = disc_left.n_components
+            disc_right = pybamm.VectorField(*[disc_right] * n)
+        else:
+            n = disc_right.n_components
+            disc_left = pybamm.VectorField(*[disc_left] * n)
+        new_comps = [
+            pybamm.simplify_if_constant(
+                symbol.create_copy(
+                    new_children=[disc_left.components[k], disc_right.components[k]]
+                )
+            )
+            for k in range(n)
+        ]
+        result = pybamm.VectorField(*new_comps)
+        result._disc_state_vector = (
+            disc_left._disc_state_vector or disc_right._disc_state_vector
+        )
+        return result
+
     def _process_symbol(self, symbol):
         """See :meth:`Discretisation.process_symbol()`."""
 
@@ -938,31 +969,9 @@ def _process_symbol(self, symbol):
                 if isinstance(disc_left, pybamm.VectorField) or isinstance(
                     disc_right, pybamm.VectorField
                 ):
-                    if isinstance(disc_left, pybamm.VectorField):
-                        n = disc_left.n_components
-                    else:
-                        n = disc_right.n_components
-                    if not isinstance(disc_right, pybamm.VectorField):
-                        disc_right = pybamm.VectorField(*[disc_right] * n)
-                    if not isinstance(disc_left, pybamm.VectorField):
-                        disc_left = pybamm.VectorField(*[disc_left] * n)
-                    new_comps = [
-                        pybamm.simplify_if_constant(
-                            symbol.create_copy(
-                                new_children=[
-                                    disc_left._components[k],
-                                    disc_right._components[k],
-                                ]
-                            )
-                        )
-                        for k in range(n)
-                    ]
-                    result = pybamm.VectorField(*new_comps)
-                    for src in (disc_left, disc_right):
-                        if hasattr(src, "_disc_state_vector"):
-                            result._disc_state_vector = src._disc_state_vector
-                            break
-                    return result
+                    return self._process_vector_field_binary(
+                        symbol, disc_left, disc_right
+                    )
 
                 return pybamm.simplify_if_constant(
                     symbol.create_copy(new_children=[disc_left, disc_right])
@@ -1106,16 +1115,16 @@ def _process_symbol(self, symbol):
                 return disc_child
             elif isinstance(symbol, pybamm.Component):
                 if not isinstance(disc_child, pybamm.VectorField):
-                    raise ValueError("Component can only be applied to a VectorField")
-                return disc_child._components[symbol.index]
+                    raise pybamm.DiscretisationError(
+                        "Component can only be applied to a VectorField"
+                    )
+                return disc_child.components[symbol.index]
             elif isinstance(symbol, pybamm.Norm):
                 if not isinstance(disc_child, pybamm.VectorField):
-                    raise ValueError("Norm can only be applied to a VectorField")
-                result = None
-                for comp in disc_child._components:
-                    sq = comp**2
-                    result = sq if result is None else result + sq
-                return result**0.5
+                    raise pybamm.DiscretisationError(
+                        "Norm can only be applied to a VectorField"
+                    )
+                return sum(c**2 for c in disc_child.components) ** 0.5
             elif isinstance(symbol, pybamm.Magnitude):
                 if not isinstance(disc_child, pybamm.VectorField):
                     raise ValueError("Magnitude can only be applied to a vector field")
@@ -1130,11 +1139,10 @@ def _process_symbol(self, symbol):
                 if isinstance(disc_child, pybamm.VectorField):
                     new_comps = [
                         symbol.create_copy(new_children=[c])
-                        for c in disc_child._components
+                        for c in disc_child.components
                     ]
                     result = pybamm.VectorField(*new_comps)
-                    if hasattr(disc_child, "_disc_state_vector"):
-                        result._disc_state_vector = disc_child._disc_state_vector
+                    result._disc_state_vector = disc_child._disc_state_vector
                     return result
                 else:
                     return symbol.create_copy(new_children=[disc_child])
@@ -1209,7 +1217,7 @@ def _process_symbol(self, symbol):
             )
 
         elif isinstance(symbol, pybamm.VectorField):
-            processed = [self.process_symbol(c) for c in symbol._components]
+            processed = [self.process_symbol(c) for c in symbol.components]
             return symbol.create_copy(new_children=processed)
 
         elif isinstance(symbol, pybamm.TensorField):
diff --git a/packages/pybamm/src/pybamm/expression_tree/unary_operators.py b/packages/pybamm/src/pybamm/expression_tree/unary_operators.py
index be3f2df585..b6a50821ba 100644
--- a/packages/pybamm/src/pybamm/expression_tree/unary_operators.py
+++ b/packages/pybamm/src/pybamm/expression_tree/unary_operators.py
@@ -1811,16 +1811,6 @@ def sign(symbol):
     return pybamm.simplify_if_constant(Sign(symbol))
 
 
-def component(symbol, index):
-    """Convenience function for creating a :class:`Component`."""
-    return Component(symbol, index)
-
-
-def norm(symbol):
-    """Convenience function for creating a :class:`Norm`."""
-    return Norm(symbol)
-
-
 def smooth_absolute_value(symbol, k):
     """
     Smooth approximation to the absolute value function. k is the smoothing parameter,
diff --git a/packages/pybamm/src/pybamm/expression_tree/vector_field.py b/packages/pybamm/src/pybamm/expression_tree/vector_field.py
index 4004c4d47c..79eddb9ccf 100644
--- a/packages/pybamm/src/pybamm/expression_tree/vector_field.py
+++ b/packages/pybamm/src/pybamm/expression_tree/vector_field.py
@@ -16,8 +16,7 @@ class VectorField(TensorField):
 
     VectorField is a convenience subclass of TensorField for rank-1 tensors
     with N >= 2 components.  Components are stored by integer index; the
-    properties ``lr_field``, ``tb_field``, and ``fb_field`` are backward-
-    compatible aliases for ``[0]``, ``[1]``, and ``[2]``.
+    properties ``lr_field`` and ``tb_field`` are aliases for ``[0]`` and ``[1]``.
 
     Parameters
     ----------
@@ -25,6 +24,9 @@ class VectorField(TensorField):
         Two or more component symbols, all sharing the same domain.
     """
 
+    # Set by discretisation for unstructured FV edge-averaging; None until then.
+    _disc_state_vector = None
+
     def __init__(self, *components):
         if len(components) < 2:
             raise ValueError(
@@ -48,29 +50,19 @@ def _from_json(cls, snippet):
     @property
     def n_components(self):
         """Number of vector components."""
-        return len(self._components)
+        return len(self.components)
 
-    # ---- backward-compatible aliases for structured-grid directions ----
+    # ---- aliases for structured-grid directions ----
 
     @property
     def lr_field(self):
         """Component 0 (left-right / x)."""
-        return self._components[0]
+        return self.components[0]
 
     @property
     def tb_field(self):
         """Component 1 (top-bottom / y)."""
-        return self._components[1]
-
-    @property
-    def fb_field(self):
-        """Component 2 (front-back / z).  Only valid for 3-component fields."""
-        if len(self._components) < 3:
-            raise AttributeError(
-                "fb_field requires at least 3 components; this VectorField has "
-                f"{len(self._components)}"
-            )
-        return self._components[2]
+        return self.components[1]
 
     def create_copy(
         self,
@@ -80,21 +72,18 @@ def create_copy(
         if new_children is None:
             new_children = [
                 c.create_copy(perform_simplifications=perform_simplifications)
-                for c in self._components
+                for c in self.components
             ]
         return VectorField(*new_children)
 
     def _to_casadi(self, t, y, y_dot, inputs, casadi_symbols):
         """See :meth:`pybamm.Symbol._to_casadi()`."""
         return casadi.vertcat(
-            *[
-                c._to_casadi_inner(t, y, y_dot, inputs, casadi_symbols)
-                for c in self._components
-            ]
+            *self._children_to_casadi(t, y, y_dot, inputs, casadi_symbols)
         )
 
     def evaluates_on_edges(self, dimension: str) -> bool:
-        statuses = [c.evaluates_on_edges(dimension) for c in self._components]
+        statuses = [c.evaluates_on_edges(dimension) for c in self.components]
         if all(statuses):
             return True
         if not any(statuses):
diff --git a/packages/pybamm/src/pybamm/solvers/processed_variable.py b/packages/pybamm/src/pybamm/solvers/processed_variable.py
index 911c01712d..263becd5ed 100644
--- a/packages/pybamm/src/pybamm/solvers/processed_variable.py
+++ b/packages/pybamm/src/pybamm/solvers/processed_variable.py
@@ -1387,6 +1387,12 @@ def process_variable(name: str, base_variables, *args, **kwargs):
         return ProcessedVariable3DSciKitFEM(name, base_variables, *args, **kwargs)
 
     if mesh and hasattr(mesh, "edges_lr") and hasattr(mesh, "edges_tb"):
+        if isinstance(base_variables[0], pybamm.VectorField):
+            raise NotImplementedError(
+                "Reading VectorField variables from a Solution is not supported "
+                "on structured 2D finite-volume meshes. Use pybamm.Component to "
+                "extract a scalar component first."
+            )
         return ProcessedVariable2DFVM(name, base_variables, *args, **kwargs)
 
     # check variable shape
diff --git a/packages/pybamm/src/pybamm/solvers/solution.py b/packages/pybamm/src/pybamm/solvers/solution.py
index 4e43b0390b..b67ce33c7e 100644
--- a/packages/pybamm/src/pybamm/solvers/solution.py
+++ b/packages/pybamm/src/pybamm/solvers/solution.py
@@ -750,7 +750,7 @@ def _update_variable(self, name: str):
                 )
             if isinstance(_var_pybamm, pybamm.VectorField):
                 comp_casadi = []
-                for k, comp in enumerate(_var_pybamm._components):
+                for k, comp in enumerate(_var_pybamm.components):
                     cc, _, _ = self._update_model_variable(
                         model,
                         comp,
diff --git a/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py b/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py
index 696f36abf0..c8f07bb9a5 100644
--- a/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py
+++ b/packages/pybamm/tests/unit/test_expression_tree/test_vector_field_and_magnitude.py
@@ -70,19 +70,33 @@ def test_component_and_norm_discretisation(self, mesh_2d):
         disc = pybamm.Discretisation(mesh_2d, spatial_methods)
         vector_field = pybamm.VectorField(pybamm.Scalar(3), pybamm.Scalar(4))
 
-        comp_0 = disc.process_symbol(pybamm.component(vector_field, 0))
-        comp_1 = disc.process_symbol(pybamm.component(vector_field, 1))
+        comp_0 = disc.process_symbol(pybamm.Component(vector_field, 0))
+        comp_1 = disc.process_symbol(pybamm.Component(vector_field, 1))
         assert comp_0.evaluate() == 3
         assert comp_1.evaluate() == 4
 
-        norm = disc.process_symbol(pybamm.norm(vector_field))
+        norm = disc.process_symbol(pybamm.Norm(vector_field))
         assert norm.evaluate() == pytest.approx(5.0)
 
-        with pytest.raises(ValueError, match=r"Component can only be applied"):
-            disc.process_symbol(pybamm.component(pybamm.Scalar(1), 0))
+        with pytest.raises(
+            pybamm.DiscretisationError, match=r"Component can only be applied"
+        ):
+            disc.process_symbol(pybamm.Component(pybamm.Scalar(1), 0))
 
-        with pytest.raises(ValueError, match=r"Norm can only be applied"):
-            disc.process_symbol(pybamm.norm(pybamm.Scalar(1)))
+        with pytest.raises(
+            pybamm.DiscretisationError, match=r"Norm can only be applied"
+        ):
+            disc.process_symbol(pybamm.Norm(pybamm.Scalar(1)))
+
+    def test_mismatched_vector_field_components(self, mesh_2d):
+        spatial_methods = {"macroscale": pybamm.FiniteVolume2D()}
+        disc = pybamm.Discretisation(mesh_2d, spatial_methods)
+        vf2 = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
+        vf3 = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2), pybamm.Scalar(3))
+        with pytest.raises(
+            pybamm.DiscretisationError, match=r"Cannot combine VectorFields"
+        ):
+            disc.process_symbol(vf2 + vf3)
 
     def test_disc_state_vector_propagation(self, mesh_2d):
         # binary and unary operators on a discretised VectorField must carry
diff --git a/packages/pybamm/tests/unit/test_solvers/test_solution.py b/packages/pybamm/tests/unit/test_solvers/test_solution.py
index 367485c265..b81cb19afe 100644
--- a/packages/pybamm/tests/unit/test_solvers/test_solution.py
+++ b/packages/pybamm/tests/unit/test_solvers/test_solution.py
@@ -513,9 +513,7 @@ def test_observable_computed_lazily(self):
         assert sol_sum._observable == expected  # and cached
 
     def test_vector_field_variable(self, mesh_2d):
-        # a VectorField model variable is processed component-by-component,
-        # storing a list of casadi functions for the downstream processed
-        # variable classes to consume
+        # structured 2D FV does not support reading VectorField solution variables
         model = pybamm.BaseModel()
         var = pybamm.Variable(
             "var", domain=["negative electrode", "separator", "positive electrode"]
@@ -528,11 +526,10 @@ def test_vector_field_variable(self, mesh_2d):
         disc.process_model(model)
         solution = pybamm.IDAKLUSolver().solve(model, np.linspace(0, 1, 5))
 
-        flux = solution["flux"]
-        assert isinstance(flux.base_variables[0], pybamm.VectorField)
-        casadi_components = flux.base_variables_casadi[0]
-        assert isinstance(casadi_components, list)
-        assert len(casadi_components) == 2
+        with pytest.raises(
+            NotImplementedError, match=r"structured 2D finite-volume meshes"
+        ):
+            solution["flux"]
 
     def test_add_solutions_different_models(self):
         # Set up first solution
diff --git a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
index 604ea629a5..151249b98c 100644
--- a/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
+++ b/packages/pybamm/tests/unit/test_spatial_methods/test_finite_volume_2d/test_tensor_field.py
@@ -201,19 +201,13 @@ def test_vectorfield_requires_two_components(self):
         with pytest.raises(ValueError, match="requires at least 2 components"):
             pybamm.VectorField(pybamm.Scalar(1))
 
-    def test_vectorfield_fb_field_three_components(self):
-        """fb_field returns 3rd component for 3-component VectorField."""
+    def test_vectorfield_three_components(self):
+        """3-component VectorField exposes components by index."""
         a, b, c = pybamm.Scalar(1), pybamm.Scalar(2), pybamm.Scalar(3)
         vf = pybamm.VectorField(a, b, c)
-        assert vf.fb_field == c
+        assert vf[2] == c
         assert vf.n_components == 3
 
-    def test_vectorfield_fb_field_raises_when_missing(self):
-        """fb_field on 2-component VectorField raises AttributeError."""
-        vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
-        with pytest.raises(AttributeError, match="fb_field requires at least 3"):
-            _ = vf.fb_field
-
     def test_vectorfield_evaluates_on_edges_all_true(self):
         """VectorField evaluates_on_edges returns True when all on edges."""
         vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
@@ -237,10 +231,10 @@ def test_vectorfield_to_casadi(self):
 class TestComponentAndNorm:
     """Tests for the Component and Norm operators on VectorFields."""
 
-    def test_component_convenience_function_and_copy(self):
-        """pybamm.component creates a Component; copying preserves the index."""
+    def test_component_copy(self):
+        """Copying a Component preserves the index."""
         vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
-        comp = pybamm.component(vf, 1)
+        comp = pybamm.Component(vf, 1)
         assert isinstance(comp, pybamm.Component)
         assert comp.index == 1
         copy = comp.create_copy()
@@ -248,10 +242,10 @@ def test_component_convenience_function_and_copy(self):
         assert copy.index == 1
         assert copy == comp
 
-    def test_norm_convenience_function_and_copy(self):
-        """pybamm.norm creates a Norm; copying preserves structure."""
+    def test_norm_copy(self):
+        """Copying a Norm preserves structure."""
         vf = pybamm.VectorField(pybamm.Scalar(1), pybamm.Scalar(2))
-        norm = pybamm.norm(vf)
+        norm = pybamm.Norm(vf)
         assert isinstance(norm, pybamm.Norm)
         copy = norm.create_copy()
         assert isinstance(copy, pybamm.Norm)

From 69b097e6b2ddf7cd0557f8bae6db59c4c83d676e Mon Sep 17 00:00:00 2001
From: Alexander Bills <alexanderallenbills@gmail.com>
Date: Mon, 3 Aug 2026 12:08:06 -0700
Subject: [PATCH 5/5] Fix VectorField _disc_state_vector propagation for
 Symbols.

`or` bool-evaluates StateVector and raises; use an explicit None check.

Co-authored-by: Cursor <cursoragent@cursor.com>
---
 .../pybamm/src/pybamm/discretisations/discretisation.py    | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/packages/pybamm/src/pybamm/discretisations/discretisation.py b/packages/pybamm/src/pybamm/discretisations/discretisation.py
index 046170af6e..cdb404505d 100644
--- a/packages/pybamm/src/pybamm/discretisations/discretisation.py
+++ b/packages/pybamm/src/pybamm/discretisations/discretisation.py
@@ -930,9 +930,10 @@ def _process_vector_field_binary(self, symbol, disc_left, disc_right):
             for k in range(n)
         ]
         result = pybamm.VectorField(*new_comps)
-        result._disc_state_vector = (
-            disc_left._disc_state_vector or disc_right._disc_state_vector
-        )
+        if disc_left._disc_state_vector is not None:
+            result._disc_state_vector = disc_left._disc_state_vector
+        else:
+            result._disc_state_vector = disc_right._disc_state_vector
         return result
 
     def _process_symbol(self, symbol):