From 17189d5c12b8d9292faaf1271d5014b5f64d8378 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Tue, 8 Nov 2022 23:20:31 -0500
Subject: [PATCH 01/12] allow concatenations of variables to be simplified

---
 pybamm/expression_tree/binary_operators.py           | 12 ++----------
 .../test_expression_tree/test_binary_operators.py    |  2 +-
 2 files changed, 3 insertions(+), 11 deletions(-)

diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 6f7a131ff4..d2e557e401 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -711,16 +711,8 @@ def _simplified_binary_broadcast_concatenation(left, right, operator):
             return left._concatenation_new_copy(
                 [operator(child, right) for child in left.orphans]
             )
-        elif (
-            isinstance(right, pybamm.Concatenation)
-            and not any(
-                isinstance(child, (pybamm.Variable, pybamm.StateVector))
-                for child in right.children
-            )
-            and (
-                all(child.is_constant() for child in left.children)
-                or all(child.is_constant() for child in right.children)
-            )
+        elif isinstance(right, pybamm.Concatenation) and not isinstance(
+            right, pybamm.ConcatenationVariable
         ):
             return left._concatenation_new_copy(
                 [
diff --git a/tests/unit/test_expression_tree/test_binary_operators.py b/tests/unit/test_expression_tree/test_binary_operators.py
index 1227ee33e0..913a943e88 100644
--- a/tests/unit/test_expression_tree/test_binary_operators.py
+++ b/tests/unit/test_expression_tree/test_binary_operators.py
@@ -557,7 +557,7 @@ def conc_broad(x, y, z):
         self.assertEqual((a + b), conc_broad(12, 14, 16))
         self.assertIsInstance((a + c), pybamm.Concatenation)
 
-        # No simplifications if are Variable or StateVector objects
+        # No simplifications if all are Variable or StateVector objects
         v = pybamm.concatenation(
             pybamm.Variable("x", "negative electrode"),
             pybamm.Variable("y", "separator"),

From 1253569122cc10e86ea2a2af005dd100848823e7 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 01:00:42 -0500
Subject: [PATCH 02/12] more simplifications again

---
 examples/scripts/DFN.py                       |  2 +-
 pybamm/discretisations/discretisation.py      |  5 ++-
 pybamm/expression_tree/averages.py            | 42 +++++++------------
 pybamm/expression_tree/binary_operators.py    |  9 ++++
 pybamm/expression_tree/broadcasts.py          |  9 +++-
 pybamm/expression_tree/unary_operators.py     | 18 +++++++-
 .../test_unary_operators.py                   | 10 +----
 7 files changed, 54 insertions(+), 41 deletions(-)

diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py
index b290b96f40..628dfbbda6 100644
--- a/examples/scripts/DFN.py
+++ b/examples/scripts/DFN.py
@@ -8,7 +8,7 @@
 pybamm.set_logging_level("INFO")
 
 # load model
-model = pybamm.lithium_ion.DFN()
+model = pybamm.lithium_ion.DFN({"working electrode": "positive", "thermal": "lumped"})
 # create geometry
 geometry = model.default_geometry
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 2068295c2b..e9d0f3f6b5 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -750,7 +750,10 @@ def process_dict(self, var_eqn_dict):
         for eqn_key, eqn in var_eqn_dict.items():
             # Broadcast if the equation evaluates to a number (e.g. Scalar)
             if np.prod(eqn.shape_for_testing) == 1 and not isinstance(eqn_key, str):
-                eqn = pybamm.FullBroadcast(eqn, broadcast_domains=eqn_key.domains)
+                if eqn_key.domain == []:
+                    eqn = eqn * pybamm.Vector([1])
+                else:
+                    eqn = pybamm.FullBroadcast(eqn, broadcast_domains=eqn_key.domains)
 
             pybamm.logger.debug("Discretise {!r}".format(eqn_key))
 
diff --git a/pybamm/expression_tree/averages.py b/pybamm/expression_tree/averages.py
index 4b95bb95ca..629be58a14 100644
--- a/pybamm/expression_tree/averages.py
+++ b/pybamm/expression_tree/averages.py
@@ -144,36 +144,22 @@ def x_average(symbol):
         else:  # pragma: no cover
             # It should be impossible to get here
             raise NotImplementedError
-    # If symbol is a concatenation of Broadcasts, its average value is the
-    # thickness-weighted average of the symbols being broadcasted
-    elif isinstance(symbol, pybamm.Concatenation) and all(
-        isinstance(child, pybamm.Broadcast) for child in symbol.children
+    # If symbol is a concatenation, its average value is the
+    # thickness-weighted average of the average of its children
+    elif isinstance(symbol, pybamm.Concatenation) and not isinstance(
+        symbol, pybamm.ConcatenationVariable
     ):
         geo = pybamm.geometric_parameters
-        l_n = geo.n.l
-        l_s = geo.s.l
-        l_p = geo.p.l
-        if symbol.domain == ["negative electrode", "separator", "positive electrode"]:
-            a, b, c = [orp.orphans[0] for orp in symbol.orphans]
-            out = (l_n * a + l_s * b + l_p * c) / (l_n + l_s + l_p)
-        elif symbol.domain == ["separator", "positive electrode"]:
-            b, c = [orp.orphans[0] for orp in symbol.orphans]
-            out = (l_s * b + l_p * c) / (l_s + l_p)
-        # To respect domains we may need to broadcast the child back out
-        child = symbol.children[0]
-        # If symbol being returned doesn't have empty domain, return it
-        if out.domain != []:
-            return out
-        # Otherwise we may need to broadcast it
-        elif child.domains["secondary"] == []:
-            return out
-        else:
-            domain = child.domains["secondary"]
-            if child.domains["tertiary"] == []:
-                return pybamm.PrimaryBroadcast(out, domain)
-            else:
-                auxiliary_domains = {"secondary": child.domains["tertiary"]}
-                return pybamm.FullBroadcast(out, domain, auxiliary_domains)
+        ls = {
+            ("negative electrode",): geo.n.l,
+            ("separator",): geo.s.l,
+            ("positive electrode",): geo.p.l,
+            ("separator", "positive electrode"): geo.s.l + geo.p.l,
+        }
+        out = sum(
+            ls[tuple(orp.domain)] * x_average(orp) for orp in symbol.orphans
+        ) / sum(ls[tuple(orp.domain)] for orp in symbol.orphans)
+        return out
     # Average of a sum is sum of averages
     elif isinstance(symbol, (pybamm.Addition, pybamm.Subtraction)):
         return _sum_of_averages(symbol, x_average)
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index d2e557e401..5712887773 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -1049,6 +1049,15 @@ def simplified_multiplication(left, right):
                     elif isinstance(right, Subtraction):
                         return (left * r_left) - (left * r_right)
 
+    # # Move constants in multiplications to the left
+    # if isinstance(left, Multiplication) and left.left.is_constant():
+    #     l_left, l_right = left.orphans
+    #     print(l_left, l_right, right)
+    #     return l_left * (l_right * right)
+    # elif isinstance(right, Multiplication) and right.left.is_constant():
+    #     r_left, r_right = right.orphans
+    #     return r_left * (left * r_right)
+
     # Cancelling out common terms
     if isinstance(left, Division):
         # Simplify (a / b) * b to a
diff --git a/pybamm/expression_tree/broadcasts.py b/pybamm/expression_tree/broadcasts.py
index 2c235ecd6d..6bfa901a42 100644
--- a/pybamm/expression_tree/broadcasts.py
+++ b/pybamm/expression_tree/broadcasts.py
@@ -87,6 +87,8 @@ def check_and_set_domains(self, child, broadcast_domain):
         # Can only do primary broadcast from current collector to electrode,
         # particle-size or particle or from electrode to particle-size or particle.
         # Note e.g. current collector to particle *is* allowed
+        if broadcast_domain == []:
+            raise pybamm.DomainError("Cannot Broadcast an object into empty domain.")
         if child.domain == []:
             pass
         elif child.domain == ["current collector"] and not (
@@ -430,7 +432,10 @@ def __init__(
 
     def check_and_set_domains(self, child, broadcast_domains):
         """See :meth:`Broadcast.check_and_set_domains`"""
-
+        if broadcast_domains["primary"] == []:
+            raise pybamm.DomainError(
+                """Cannot do full broadcast to an empty primary domain"""
+            )
         # Variables on the current collector can only be broadcast to 'primary'
         if child.domain == ["current collector"]:
             raise pybamm.DomainError(
@@ -544,6 +549,8 @@ def full_like(symbols, fill_value):
         return array_type(entries, domains=sum_symbol.domains)
 
     except NotImplementedError:
+        if sum_symbol.shape_for_testing == (1, 1):
+            return pybamm.Scalar(fill_value)
         if sum_symbol.evaluates_on_edges("primary"):
             return FullBroadcastToEdges(
                 fill_value, broadcast_domains=sum_symbol.domains
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index c6ab47c116..5b4ce7af06 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -1084,7 +1084,10 @@ def grad(symbol):
     """
     # Gradient of a broadcast is zero
     if isinstance(symbol, pybamm.PrimaryBroadcast):
-        new_child = pybamm.PrimaryBroadcast(0, symbol.child.domain)
+        if symbol.child.domain == []:
+            new_child = pybamm.Scalar(0)
+        else:
+            new_child = pybamm.PrimaryBroadcast(0, symbol.child.domain)
         return pybamm.PrimaryBroadcastToEdges(new_child, symbol.domain)
     elif isinstance(symbol, pybamm.FullBroadcast):
         return pybamm.FullBroadcastToEdges(0, broadcast_domains=symbol.domains)
@@ -1110,7 +1113,10 @@ def div(symbol):
     """
     # Divergence of a broadcast is zero
     if isinstance(symbol, pybamm.PrimaryBroadcastToEdges):
-        new_child = pybamm.PrimaryBroadcast(0, symbol.child.domain)
+        if symbol.child.domain == []:
+            new_child = pybamm.Scalar(0)
+        else:
+            new_child = pybamm.PrimaryBroadcast(0, symbol.child.domain)
         return pybamm.PrimaryBroadcast(new_child, symbol.domain)
     # Divergence commutes with Negate operator
     if isinstance(symbol, pybamm.Negate):
@@ -1245,6 +1251,14 @@ def boundary_value(symbol, side):
 
 def sign(symbol):
     """Returns a :class:`Sign` object."""
+    if isinstance(symbol, pybamm.Broadcast):
+        # Move sign inside the broadcast
+        # Apply recursively
+        return symbol._unary_new_copy(sign(symbol.orphans[0]))
+    elif isinstance(symbol, pybamm.Concatenation) and not isinstance(
+        symbol, pybamm.ConcatenationVariable
+    ):
+        return pybamm.concatenation(*[sign(child) for child in symbol.orphans])
     return pybamm.simplify_if_constant(Sign(symbol))
 
 
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index dd587228db..37937fbead 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -147,10 +147,7 @@ def test_gradient(self):
         # gradient of broadcast should return broadcasted zero
         a = pybamm.PrimaryBroadcast(pybamm.Variable("a"), "test domain")
         grad = pybamm.grad(a)
-        self.assertIsInstance(grad, pybamm.PrimaryBroadcastToEdges)
-        self.assertIsInstance(grad.child, pybamm.PrimaryBroadcast)
-        self.assertIsInstance(grad.child.child, pybamm.Scalar)
-        self.assertEqual(grad.child.child.value, 0)
+        self.assertEqual(grad, pybamm.PrimaryBroadcastToEdges(0, "test domain"))
 
         # otherwise gradient should work
         a = pybamm.Symbol("a", domain="test domain")
@@ -175,10 +172,7 @@ def test_div(self):
         # divergence of broadcast should return broadcasted zero
         a = pybamm.PrimaryBroadcastToEdges(pybamm.Variable("a"), "test domain")
         div = pybamm.div(a)
-        self.assertIsInstance(div, pybamm.PrimaryBroadcast)
-        self.assertIsInstance(div.child, pybamm.PrimaryBroadcast)
-        self.assertIsInstance(div.child.child, pybamm.Scalar)
-        self.assertEqual(div.child.child.value, 0)
+        self.assertEqual(div, pybamm.PrimaryBroadcast(0, "test domain"))
 
         # otherwise divergence should work
         a = pybamm.Symbol("a", domain="test domain")

From e388a89eb28219a0d59e5d4433446a26349ab857 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 10:25:53 -0500
Subject: [PATCH 03/12] coverage

---
 examples/scripts/compare_lithium_ion.py       |  7 +++---
 pybamm/discretisations/discretisation.py      | 10 +++------
 pybamm/expression_tree/symbol.py              |  2 +-
 pybamm/expression_tree/unary_operators.py     |  3 +++
 .../test_expression_tree/test_broadcasts.py   |  9 ++++++++
 .../test_unary_operators.py                   | 22 +++++++++++++++++++
 6 files changed, 42 insertions(+), 11 deletions(-)

diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py
index 6108036b9b..4538f8b2f1 100644
--- a/examples/scripts/compare_lithium_ion.py
+++ b/examples/scripts/compare_lithium_ion.py
@@ -7,10 +7,10 @@
 
 # load models
 models = [
-    pybamm.lithium_ion.SPM(),
+    # pybamm.lithium_ion.SPM(),
     pybamm.lithium_ion.SPMe(),
-    pybamm.lithium_ion.DFN(),
-    pybamm.lithium_ion.NewmanTobias(),
+    # pybamm.lithium_ion.DFN(),
+    # pybamm.lithium_ion.NewmanTobias(),
 ]
 
 # create and run simulations
@@ -18,6 +18,7 @@
 for model in models:
     sim = pybamm.Simulation(model)
     sim.solve([0, 3600])
+    sim.built_model.variables["Terminal voltage [V]"].render()
     sims.append(sim)
 
 # plot
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index e9d0f3f6b5..218b5be678 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -900,13 +900,9 @@ def _process_symbol(self, symbol):
             elif isinstance(symbol, pybamm.Broadcast):
                 # Broadcast new_child to the domain specified by symbol.domain
                 # Different discretisations may broadcast differently
-                if symbol.domain == []:
-                    out = disc_child * pybamm.Vector([1])
-                else:
-                    out = spatial_method.broadcast(
-                        disc_child, symbol.domains, symbol.broadcast_type
-                    )
-                return out
+                return spatial_method.broadcast(
+                    disc_child, symbol.domains, symbol.broadcast_type
+                )
 
             elif isinstance(symbol, pybamm.DeltaFunction):
                 return spatial_method.delta_function(symbol, disc_child)
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 5ec7dbf1f4..ddd959ab21 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -615,7 +615,7 @@ def __abs__(self):
         elif isinstance(self, pybamm.Broadcast):
             # Move absolute value inside the broadcast
             # Apply recursively
-            abs_self_not_broad = pybamm.simplify_if_constant(abs(self.orphans[0]))
+            abs_self_not_broad = abs(self.orphans[0])
             return self._unary_new_copy(abs_self_not_broad)
         else:
             k = pybamm.settings.abs_smoothing
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index 5b4ce7af06..860fb1dbaa 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -965,6 +965,9 @@ def __init__(self, children, initial_condition):
     def _unary_new_copy(self, child):
         return self.__class__(child, self.initial_condition)
 
+    def is_constant(self):
+        return False
+
 
 class BoundaryGradient(BoundaryOperator):
     """
diff --git a/tests/unit/test_expression_tree/test_broadcasts.py b/tests/unit/test_expression_tree/test_broadcasts.py
index 300deb9710..f9500a6f90 100644
--- a/tests/unit/test_expression_tree/test_broadcasts.py
+++ b/tests/unit/test_expression_tree/test_broadcasts.py
@@ -52,6 +52,10 @@ def test_primary_broadcast(self):
         )
 
         a = pybamm.Symbol("a", domain="current collector")
+        with self.assertRaisesRegex(
+            pybamm.DomainError, "Cannot Broadcast an object into empty domain"
+        ):
+            pybamm.PrimaryBroadcast(a, [])
         with self.assertRaisesRegex(
             pybamm.DomainError, "Primary broadcast from current collector"
         ):
@@ -209,6 +213,11 @@ def test_full_broadcast(self):
             ),
         )
 
+        with self.assertRaisesRegex(
+            pybamm.DomainError, "Cannot do full broadcast to an empty primary domain"
+        ):
+            pybamm.FullBroadcast(a, [])
+
     def test_full_broadcast_number(self):
         broad_a = pybamm.FullBroadcast(1, ["negative electrode"], None)
         self.assertEqual(broad_a.name, "broadcast")
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 37937fbead..9dd36d4d7d 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -103,6 +103,18 @@ def test_sign(self):
             np.diag(signb.evaluate().toarray()), [-1, -1, 0, 1, 1]
         )
 
+        broad = pybamm.PrimaryBroadcast(-4, "test domain")
+        self.assertEqual(pybamm.sign(broad), pybamm.PrimaryBroadcast(-1, "test domain"))
+
+        conc = pybamm.Concatenation(broad, pybamm.PrimaryBroadcast(2, "another domain"))
+        self.assertEqual(
+            pybamm.sign(conc),
+            pybamm.Concatenation(
+                pybamm.PrimaryBroadcast(-1, "test domain"),
+                pybamm.PrimaryBroadcast(1, "another domain"),
+            ),
+        )
+
     def test_floor(self):
         a = pybamm.Symbol("a")
         floora = pybamm.Floor(a)
@@ -173,6 +185,16 @@ def test_div(self):
         a = pybamm.PrimaryBroadcastToEdges(pybamm.Variable("a"), "test domain")
         div = pybamm.div(a)
         self.assertEqual(div, pybamm.PrimaryBroadcast(0, "test domain"))
+        a = pybamm.PrimaryBroadcastToEdges(
+            pybamm.Variable("a", "some domain"), "test domain"
+        )
+        div = pybamm.div(a)
+        self.assertEqual(
+            div,
+            pybamm.PrimaryBroadcast(
+                pybamm.PrimaryBroadcast(0, "some domain"), "test domain"
+            ),
+        )
 
         # otherwise divergence should work
         a = pybamm.Symbol("a", domain="test domain")

From 898d32f203fae300f290b6ed64c511a585267a72 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 11:13:05 -0500
Subject: [PATCH 04/12] make electrolyte concentration cleaner

---
 examples/scripts/DFN.py                       |  2 +-
 examples/scripts/compare_lithium_ion.py       |  7 ++---
 .../base_electrolyte_diffusion.py             | 29 +++++++++++++++----
 .../electrolyte_diffusion/full_diffusion.py   | 12 ++++++--
 4 files changed, 38 insertions(+), 12 deletions(-)

diff --git a/examples/scripts/DFN.py b/examples/scripts/DFN.py
index 628dfbbda6..b290b96f40 100644
--- a/examples/scripts/DFN.py
+++ b/examples/scripts/DFN.py
@@ -8,7 +8,7 @@
 pybamm.set_logging_level("INFO")
 
 # load model
-model = pybamm.lithium_ion.DFN({"working electrode": "positive", "thermal": "lumped"})
+model = pybamm.lithium_ion.DFN()
 # create geometry
 geometry = model.default_geometry
 
diff --git a/examples/scripts/compare_lithium_ion.py b/examples/scripts/compare_lithium_ion.py
index 4538f8b2f1..6108036b9b 100644
--- a/examples/scripts/compare_lithium_ion.py
+++ b/examples/scripts/compare_lithium_ion.py
@@ -7,10 +7,10 @@
 
 # load models
 models = [
-    # pybamm.lithium_ion.SPM(),
+    pybamm.lithium_ion.SPM(),
     pybamm.lithium_ion.SPMe(),
-    # pybamm.lithium_ion.DFN(),
-    # pybamm.lithium_ion.NewmanTobias(),
+    pybamm.lithium_ion.DFN(),
+    pybamm.lithium_ion.NewmanTobias(),
 ]
 
 # create and run simulations
@@ -18,7 +18,6 @@
 for model in models:
     sim = pybamm.Simulation(model)
     sim.solve([0, 3600])
-    sim.built_model.variables["Terminal voltage [V]"].render()
     sims.append(sim)
 
 # plot
diff --git a/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py b/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py
index 6deb282e8f..0af2ada0e0 100644
--- a/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py
+++ b/pybamm/models/submodels/electrolyte_diffusion/base_electrolyte_diffusion.py
@@ -37,16 +37,17 @@ def _get_standard_concentration_variables(self, c_e_dict):
             electrolyte.
         """
 
-        c_e_typ = self.param.c_e_typ
         c_e = pybamm.concatenation(*c_e_dict.values())
         # Override print_name
         c_e.print_name = "c_e"
 
-        variables = {
-            "Electrolyte concentration": c_e,
-            "X-averaged electrolyte concentration": pybamm.x_average(c_e),
-        }
+        variables = self._get_standard_domain_concentration_variables(c_e_dict)
+        variables.update(self._get_standard_whole_cell_concentration_variables(c_e))
+        return variables
 
+    def _get_standard_domain_concentration_variables(self, c_e_dict):
+        c_e_typ = self.param.c_e_typ
+        variables = {}
         # Case where an electrode is not included (half-cell)
         if "negative electrode" not in self.options.whole_cell_domains:
             c_e_s = c_e_dict["separator"]
@@ -75,6 +76,24 @@ def _get_standard_concentration_variables(self, c_e_dict):
 
         return variables
 
+    def _get_standard_whole_cell_concentration_variables(self, c_e):
+        c_e_typ = self.param.c_e_typ
+
+        variables = {
+            "Electrolyte concentration": c_e,
+            "X-averaged electrolyte concentration": pybamm.x_average(c_e),
+        }
+        variables_nondim = variables.copy()
+        for name, var in variables_nondim.items():
+            variables.update(
+                {
+                    f"{name} [mol.m-3]": c_e_typ * var,
+                    f"{name} [Molar]": c_e_typ * var / 1000,
+                }
+            )
+
+        return variables
+
     def _get_standard_porosity_times_concentration_variables(self, eps_c_e_dict):
         eps_c_e = pybamm.concatenation(*eps_c_e_dict.values())
         variables = {"Porosity times concentration": eps_c_e}
diff --git a/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py b/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
index a7cd2bd6aa..411cefb1fb 100644
--- a/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
+++ b/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
@@ -53,10 +53,15 @@ def get_coupled_variables(self, variables):
             c_e_k = eps_c_e_k / eps_k
             c_e_dict[domain] = c_e_k
 
-        variables.update(self._get_standard_concentration_variables(c_e_dict))
+        variables.update(self._get_standard_domain_concentration_variables(c_e_dict))
+
+        c_e = variables["Porosity times concentration"] / variables["Porosity"]
+        variables.update(self._get_standard_whole_cell_concentration_variables(c_e))
+        variables["Electrolyte concentration concatenation"] = pybamm.concatenation(
+            *c_e_dict.values()
+        )
 
         # Whole domain
-        c_e = variables["Electrolyte concentration"]
         tor = variables["Electrolyte transport efficiency"]
         i_e = variables["Electrolyte current density"]
         v_box = variables["Volume-averaged velocity"]
@@ -102,6 +107,7 @@ def set_initial_conditions(self, variables):
     def set_boundary_conditions(self, variables):
         param = self.param
         c_e = variables["Electrolyte concentration"]
+        c_e_conc = variables["Electrolyte concentration concatenation"]
         T = variables["Cell temperature"]
         tor = variables["Electrolyte transport efficiency"]
         i_boundary_cc = variables["Current collector current density"]
@@ -132,6 +138,8 @@ def flux_bc(side):
         #     # right bc at separator/cathode interface
         #     rbc = flux_bc("right")
 
+        # add boundary conditions to both forms of the concentration
         self.boundary_conditions = {
             c_e: {"left": (lbc, "Neumann"), "right": (rbc, "Neumann")},
+            c_e_conc: {"left": (lbc, "Neumann"), "right": (rbc, "Neumann")},
         }

From 796bf809d3643e4de6b41fb014d24fd25773b1d2 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 14:02:57 -0500
Subject: [PATCH 05/12] fix test and memoize submeshes

---
 .../compare_comsol/compare_comsol_DFN.py      |  4 +-
 pybamm/discretisations/discretisation.py      |  8 +-
 pybamm/expression_tree/concatenations.py      |  4 +-
 pybamm/meshes/meshes.py                       | 27 ++++--
 .../electrolyte_diffusion/full_diffusion.py   |  6 +-
 pybamm/spatial_methods/finite_volume.py       | 36 ++++----
 pybamm/spatial_methods/spatial_method.py      | 10 +--
 pybamm/spatial_methods/spectral_volume.py     | 13 ++-
 .../test_models/standard_output_tests.py      |  4 +-
 .../test_interface/test_butler_volmer.py      |  2 +-
 .../test_interface/test_lead_acid.py          |  2 +-
 .../test_interface/test_lithium_ion.py        |  2 +-
 .../test_finite_volume.py                     | 10 +--
 .../test_spectral_volume.py                   | 10 +--
 .../test_discretisation.py                    | 90 +++++++++----------
 .../test_operations/test_evaluate_julia.py    | 14 ++-
 tests/unit/test_meshes/test_meshes.py         |  3 +-
 .../test_lead_acid_parameters.py              |  6 +-
 tests/unit/test_plotting/test_quick_plot.py   |  2 +-
 .../unit/test_solvers/test_scikits_solvers.py | 12 +--
 tests/unit/test_solvers/test_scipy_solver.py  |  6 +-
 .../test_base_spatial_method.py               |  4 +-
 .../test_finite_volume/test_extrapolation.py  |  4 +-
 .../test_finite_volume/test_finite_volume.py  | 22 ++---
 .../test_ghost_nodes_and_neumann.py           |  8 +-
 .../test_grad_div_shapes.py                   | 74 ++++++++-------
 .../test_finite_volume/test_integration.py    | 90 +++++++++----------
 .../test_spectral_volume.py                   | 74 ++++++++-------
 .../test_zero_dimensional_method.py           |  4 +-
 29 files changed, 273 insertions(+), 278 deletions(-)

diff --git a/examples/scripts/compare_comsol/compare_comsol_DFN.py b/examples/scripts/compare_comsol/compare_comsol_DFN.py
index 1d17f7ba33..e5a9930a56 100644
--- a/examples/scripts/compare_comsol/compare_comsol_DFN.py
+++ b/examples/scripts/compare_comsol/compare_comsol_DFN.py
@@ -78,7 +78,7 @@ def get_interp_fun(variable_name, domain):
         comsol_x = comsol_variables["x"]
 
     # Make sure to use dimensional space
-    pybamm_x = mesh.combine_submeshes(*domain).nodes * L_x
+    pybamm_x = mesh[domain].nodes * L_x
     variable = interp.interp1d(comsol_x, variable, axis=0)(pybamm_x)
 
     fun = pybamm.Interpolant(
@@ -88,7 +88,7 @@ def get_interp_fun(variable_name, domain):
     )
 
     fun.domains = {"primary": domain}
-    fun.mesh = mesh.combine_submeshes(*domain)
+    fun.mesh = mesh.[domain]
     fun.secondary_mesh = None
     return fun
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 218b5be678..1243ee3586 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -787,14 +787,14 @@ def process_symbol(self, symbol):
 
             # Assign mesh as an attribute to the processed variable
             if symbol.domain != []:
-                discretised_symbol.mesh = self.mesh.combine_submeshes(*symbol.domain)
+                discretised_symbol.mesh = self.mesh[symbol.domain]
             else:
                 discretised_symbol.mesh = None
             # Assign secondary mesh
             if symbol.domains["secondary"] != []:
-                discretised_symbol.secondary_mesh = self.mesh.combine_submeshes(
-                    *symbol.domains["secondary"]
-                )
+                discretised_symbol.secondary_mesh = self.mesh[
+                    symbol.domains["secondary"]
+                ]
             else:
                 discretised_symbol.secondary_mesh = None
             return discretised_symbol
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index b3abef26f6..01b2ea870b 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -260,7 +260,7 @@ def _get_auxiliary_domain_repeats(self, auxiliary_domains):
         mesh_pts = 1
         for level, dom in auxiliary_domains.items():
             if level != "primary" and dom != []:
-                mesh_pts *= self.full_mesh.combine_submeshes(*dom).npts
+                mesh_pts *= self.full_mesh[dom].npts
         return mesh_pts
 
     @property
@@ -437,6 +437,8 @@ def simplified_concatenation(*children):
 def concatenation(*children):
     """Helper function to create concatenations."""
     # TODO: add option to turn off simplifications
+    if len(children) == 3 and children[-1].domain == ["current collector"]:
+        print("here")
     return simplified_concatenation(*children)
 
 
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 1a96137f9b..9f5208e2ab 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -117,6 +117,22 @@ def __init__(self, geometry, submesh_types, var_pts):
         # add ghost meshes
         self.add_ghost_meshes()
 
+    def __getitem__(self, domains):
+        if isinstance(domains, str):
+            domains = (domains,)
+        domains = tuple(domains)
+        try:
+            return super().__getitem__(domains)
+        except KeyError:
+            value = self.combine_submeshes(*domains)
+            self[domains] = value
+            return value
+
+    def __setitem__(self, domains, value):
+        if isinstance(domains, str):
+            domains = (domains,)
+        super().__setitem__(domains, value)
+
     def combine_submeshes(self, *submeshnames):
         """Combine submeshes into a new submesh, using self.submeshclass
         Raises pybamm.DomainError if submeshes to be combined do not match up (edges are
@@ -159,7 +175,6 @@ def combine_submeshes(self, *submeshnames):
         submesh.internal_boundaries = [
             self[submeshname].edges[0] for submeshname in submeshnames[1:]
         ]
-
         return submesh
 
     def add_ghost_meshes(self):
@@ -172,22 +187,24 @@ def add_ghost_meshes(self):
         submeshes = [
             (domain, submesh)
             for domain, submesh in self.items()
-            if not isinstance(submesh, (pybamm.SubMesh0D, pybamm.ScikitSubMesh2D))
+            if (
+                len(domain) == 1
+                and not isinstance(submesh, (pybamm.SubMesh0D, pybamm.ScikitSubMesh2D))
+            )
         ]
         for domain, submesh in submeshes:
-
             edges = submesh.edges
 
             # left ghost cell: two edges, one node, to the left of existing submesh
             lgs_edges = np.array([2 * edges[0] - edges[1], edges[0]])
-            self[domain + "_left ghost cell"] = pybamm.SubMesh1D(
+            self[domain[0] + "_left ghost cell"] = pybamm.SubMesh1D(
                 lgs_edges, submesh.coord_sys
             )
 
             # right ghost cell: two edges, one node, to the right of
             # existing submesh
             rgs_edges = np.array([edges[-1], 2 * edges[-1] - edges[-2]])
-            self[domain + "_right ghost cell"] = pybamm.SubMesh1D(
+            self[domain[0] + "_right ghost cell"] = pybamm.SubMesh1D(
                 rgs_edges, submesh.coord_sys
             )
 
diff --git a/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py b/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
index 411cefb1fb..4226f75689 100644
--- a/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
+++ b/pybamm/models/submodels/electrolyte_diffusion/full_diffusion.py
@@ -53,13 +53,13 @@ def get_coupled_variables(self, variables):
             c_e_k = eps_c_e_k / eps_k
             c_e_dict[domain] = c_e_k
 
+        variables["Electrolyte concentration concatenation"] = pybamm.concatenation(
+            *c_e_dict.values()
+        )
         variables.update(self._get_standard_domain_concentration_variables(c_e_dict))
 
         c_e = variables["Porosity times concentration"] / variables["Porosity"]
         variables.update(self._get_standard_whole_cell_concentration_variables(c_e))
-        variables["Electrolyte concentration concatenation"] = pybamm.concatenation(
-            *c_e_dict.values()
-        )
 
         # Whole domain
         tor = variables["Electrolyte transport efficiency"]
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index 25c61c71c8..8b069ff626 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -57,7 +57,7 @@ def spatial_variable(self, symbol):
         :class:`pybamm.Vector`
             Contains the discretised spatial variable
         """
-        symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
+        symbol_mesh = self.mesh[symbol.domain]
         repeats = self._get_auxiliary_domain_repeats(symbol.domains)
         if symbol.evaluates_on_edges("primary"):
             entries = np.tile(symbol_mesh.edges, repeats)
@@ -137,7 +137,7 @@ def gradient_matrix(self, domain, domains):
             The (sparse) finite volume gradient matrix for the domain
         """
         # Create appropriate submesh by combining submeshes in primary domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self.mesh[domain]
 
         # Create 1D matrix using submesh
         n = submesh.npts
@@ -160,7 +160,7 @@ def divergence(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the divergence operator.
         See :meth:`pybamm.SpatialMethod.divergence`
         """
-        submesh = self.mesh.combine_submeshes(*symbol.domain)
+        submesh = self.mesh[symbol.domain]
 
         divergence_matrix = self.divergence_matrix(symbol.domains)
 
@@ -195,7 +195,7 @@ def divergence_matrix(self, domains):
             The (sparse) finite volume divergence matrix for the domain
         """
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domains["primary"])
+        submesh = self.mesh[domains["primary"]]
 
         # check coordinate system
         if submesh.coord_sys in ["cylindrical polar", "spherical polar"]:
@@ -276,7 +276,7 @@ def definite_integral_matrix(
             )
 
         domain = child.domains[integration_dimension]
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self.mesh[domain]
 
         # check coordinate system
         if submesh.coord_sys in ["cylindrical polar", "spherical polar"]:
@@ -291,7 +291,7 @@ def definite_integral_matrix(
 
         if integration_dimension == "primary":
             # Create appropriate submesh by combining submeshes in domain
-            submesh = self.mesh.combine_submeshes(*domains["primary"])
+            submesh = self.mesh[domains["primary"]]
 
             # Create vector of ones for primary domain submesh
 
@@ -306,7 +306,7 @@ def definite_integral_matrix(
             matrix = kron(eye(second_dim_repeats), d_edges)
         elif integration_dimension == "secondary":
             # Create appropriate submesh by combining submeshes in domain
-            primary_submesh = self.mesh.combine_submeshes(*domains["primary"])
+            primary_submesh = self.mesh[domains["primary"]]
 
             # Create matrix which integrates in the secondary dimension
             # Different number of edges depending on whether child evaluates on edges
@@ -348,7 +348,7 @@ def indefinite_integral(self, child, discretised_child, direction):
             # the case where child evaluates on edges
             # If it becomes necessary to implement this, will need to think about what
             # the cylindrical/spherical polar indefinite integral should be
-            submesh = self.mesh.combine_submeshes(*child.domain)
+            submesh = self.mesh[child.domain]
             if submesh.coord_sys in ["cylindrical polar", "spherical polar"]:
                 raise NotImplementedError(
                     f"Indefinite integral on a {submesh.coord_sys} domain is not "
@@ -438,7 +438,7 @@ def indefinite_integral_matrix_edges(self, domains, direction):
         """
 
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domains["primary"])
+        submesh = self.mesh[domains["primary"]]
         n = submesh.npts
         second_dim_repeats = self._get_auxiliary_domain_repeats(domains)
 
@@ -488,7 +488,7 @@ def indefinite_integral_matrix_nodes(self, domains, direction):
         """
 
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domains["primary"])
+        submesh = self.mesh[domains["primary"]]
         n = submesh.npts
         second_dim_repeats = self._get_auxiliary_domain_repeats(domains)
 
@@ -518,7 +518,7 @@ def delta_function(self, symbol, discretised_symbol):
         See :meth:`pybamm.SpatialMethod.delta_function`
         """
         # Find the number of submeshes
-        submesh = self.mesh.combine_submeshes(*symbol.domain)
+        submesh = self.mesh[symbol.domain]
 
         prim_pts = submesh.npts
         second_dim_repeats = self._get_auxiliary_domain_repeats(symbol.domains)
@@ -635,7 +635,7 @@ def add_ghost_nodes(self, symbol, discretised_symbol, bcs):
         """
         # get relevant grid points
         domain = symbol.domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self.mesh[domain]
 
         # Prepare sizes and empty bcs_vector
         n = submesh.npts
@@ -757,7 +757,7 @@ def add_neumann_values(self, symbol, discretised_gradient, bcs, domain):
 
         """
         # get relevant grid points
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self.mesh[domain]
 
         # Prepare sizes and empty bcs_vector
         n = submesh.npts - 1
@@ -849,7 +849,7 @@ def boundary_value_or_flux(self, symbol, discretised_child, bcs=None):
         """
 
         # Find the number of submeshes
-        submesh = self.mesh.combine_submeshes(*discretised_child.domain)
+        submesh = self.mesh[discretised_child.domain]
 
         prim_pts = submesh.npts
         repeats = self._get_auxiliary_domain_repeats(discretised_child.domains)
@@ -1136,7 +1136,7 @@ def concatenation(self, disc_children):
         See :meth:`pybamm.SpatialMethod.concatenation`
         """
         for idx, child in enumerate(disc_children):
-            submesh = self.mesh.combine_submeshes(*child.domain)
+            submesh = self.mesh[child.domain]
             repeats = self._get_auxiliary_domain_repeats(child.domains)
             n_nodes = len(submesh.nodes) * repeats
             n_edges = len(submesh.edges) * repeats
@@ -1203,7 +1203,7 @@ def shift(self, discretised_symbol, shift_key, method):
         def arithmetic_mean(array):
             """Calculate the arithmetic mean of an array using matrix multiplication"""
             # Create appropriate submesh by combining submeshes in domain
-            submesh = self.mesh.combine_submeshes(*array.domain)
+            submesh = self.mesh[array.domain]
 
             # Create 1D matrix using submesh
             n = submesh.npts
@@ -1261,7 +1261,7 @@ def harmonic_mean(array):
             approximation to the diffusion equation." (2012).
             """
             # Create appropriate submesh by combining submeshes in domain
-            submesh = self.mesh.combine_submeshes(*array.domain)
+            submesh = self.mesh[array.domain]
 
             # Get second dimension length for use later
             second_dim_repeats = self._get_auxiliary_domain_repeats(
@@ -1398,7 +1398,7 @@ def upwind_or_downwind(self, symbol, discretised_symbol, bcs, direction):
         direction : str
             Direction in which to apply the operator (upwind or downwind)
         """
-        submesh = self.mesh.combine_submeshes(*symbol.domain)
+        submesh = self.mesh[symbol.domain]
         n = submesh.npts
 
         if symbol not in bcs:
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
index 267bdda3e8..6826414bab 100644
--- a/pybamm/spatial_methods/spatial_method.py
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -47,7 +47,7 @@ def _get_auxiliary_domain_repeats(self, domains):
         mesh_pts = 1
         for level, dom in domains.items():
             if level != "primary" and dom != []:
-                mesh_pts *= self.mesh.combine_submeshes(*dom).npts
+                mesh_pts *= self.mesh[dom].npts
         return mesh_pts
 
     @property
@@ -70,7 +70,7 @@ def spatial_variable(self, symbol):
         :class:`pybamm.Vector`
             Contains the discretised spatial variable
         """
-        symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
+        symbol_mesh = self.mesh[symbol.domain]
         repeats = self._get_auxiliary_domain_repeats(symbol.domains)
         if symbol.evaluates_on_edges("primary"):
             entries = np.tile(symbol_mesh.edges, repeats)
@@ -99,7 +99,7 @@ def broadcast(self, symbol, domains, broadcast_type):
             The discretised symbol of the correct size for the spatial method
         """
         domain = domains["primary"]
-        primary_domain_size = self.mesh.combine_submeshes(*domain).npts
+        primary_domain_size = self.mesh[domain].npts
         secondary_domain_size = self._get_auxiliary_domain_repeats(
             {"secondary": domains["secondary"]}
         )
@@ -403,8 +403,8 @@ def mass_matrix(self, symbol, boundary_conditions):
         # to account for Dirichlet boundary conditions. Here, we just have the default
         # behaviour that the mass matrix is the identity.
 
-        # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*symbol.domain)
+        # Get submesh
+        submesh = self.mesh[symbol.domain]
 
         # Get number of points in primary dimension
         n = submesh.npts
diff --git a/pybamm/spatial_methods/spectral_volume.py b/pybamm/spatial_methods/spectral_volume.py
index dd45fbaf2c..17fe70e040 100644
--- a/pybamm/spatial_methods/spectral_volume.py
+++ b/pybamm/spatial_methods/spectral_volume.py
@@ -160,8 +160,7 @@ def cv_boundary_reconstruction_matrix(self, domains):
         :class:`pybamm.Matrix`
             The (sparse) CV reconstruction matrix for the domain
         """
-        # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domains["primary"])
+        submesh = self.mesh[domains["primary"]]
 
         # Obtain the basic reconstruction matrix.
         recon_sub_matrix = self.cv_boundary_reconstruction_sub_matrix()
@@ -316,8 +315,7 @@ def gradient_matrix(self, domain, domains):
         :class:`pybamm.Matrix`
             The (sparse) Spectral Volume gradient matrix for the domain
         """
-        # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self.mesh[domain]
 
         # Obtain the Chebyshev differentiation matrix.
         # Flip it, since it is defined for the Chebyshev
@@ -400,8 +398,7 @@ def penalty_matrix(self, domains):
         :class:`pybamm.Matrix`
             The (sparse) Spectral Volume penalty matrix for the domain
         """
-        # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domains["primary"])
+        submesh = self.mesh[domains["primary"]]
 
         # Create 1D matrix using submesh
         n = submesh.npts
@@ -523,7 +520,7 @@ def replace_dirichlet_values(self, symbol, discretised_symbol, bcs):
         """
         # get relevant grid points
         domain = symbol.domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self.mesh[domain]
 
         # Prepare sizes
         n = (submesh.npts // self.order) * (self.order + 1)
@@ -617,7 +614,7 @@ def replace_neumann_values(self, symbol, discretised_gradient, bcs):
         """
         # get relevant grid points
         domain = symbol.domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self.mesh[domain]
 
         # Prepare sizes
         n = submesh.npts + 1
diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 641caec961..295837a63c 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -83,11 +83,11 @@ def __init__(self, model, param, disc, solution, operating_condition):
         self.x_s = disc.mesh["separator"].nodes * L_x
         self.x_p = disc.mesh["positive electrode"].nodes * L_x
         whole_cell = ["negative electrode", "separator", "positive electrode"]
-        self.x = disc.mesh.combine_submeshes(*whole_cell).nodes * L_x
+        self.x = disc.mesh[whole_cell].nodes * L_x
         self.x_n_edge = disc.mesh["negative electrode"].edges * L_x
         self.x_s_edge = disc.mesh["separator"].edges * L_x
         self.x_p_edge = disc.mesh["positive electrode"].edges * L_x
-        self.x_edge = disc.mesh.combine_submeshes(*whole_cell).edges * L_x
+        self.x_edge = disc.mesh[whole_cell].edges * L_x
 
         if isinstance(self.model, pybamm.lithium_ion.BaseModel):
             R_n_typ = model.length_scales["negative particle"].evaluate()
diff --git a/tests/integration/test_models/test_submodels/test_interface/test_butler_volmer.py b/tests/integration/test_models/test_submodels/test_interface/test_butler_volmer.py
index 8ee3e1472a..b2a6966348 100644
--- a/tests/integration/test_models/test_submodels/test_interface/test_butler_volmer.py
+++ b/tests/integration/test_models/test_submodels/test_interface/test_butler_volmer.py
@@ -223,7 +223,7 @@ def test_discretisation(self):
 
         # test concatenated butler-volmer
         whole_cell = ["negative electrode", "separator", "positive electrode"]
-        whole_cell_mesh = disc.mesh.combine_submeshes(*whole_cell)
+        whole_cell_mesh = disc.mesh[whole_cell]
         self.assertEqual(j.evaluate(None, y).shape, (whole_cell_mesh.npts, 1))
 
     def test_diff_c_e_lead_acid(self):
diff --git a/tests/integration/test_models/test_submodels/test_interface/test_lead_acid.py b/tests/integration/test_models/test_submodels/test_interface/test_lead_acid.py
index 54bd7c678b..de96380a15 100644
--- a/tests/integration/test_models/test_submodels/test_interface/test_lead_acid.py
+++ b/tests/integration/test_models/test_submodels/test_interface/test_lead_acid.py
@@ -83,7 +83,7 @@ def test_discretisation_main_reaction(self):
 
         # Test
         whole_cell = ["negative electrode", "separator", "positive electrode"]
-        submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
         y = submesh.nodes**2
         # should evaluate to vectors with the right shape
         self.assertEqual(j0_n.evaluate(y=y).shape, (mesh["negative electrode"].npts, 1))
diff --git a/tests/integration/test_models/test_submodels/test_interface/test_lithium_ion.py b/tests/integration/test_models/test_submodels/test_interface/test_lithium_ion.py
index fff744c42b..c1970d527e 100644
--- a/tests/integration/test_models/test_submodels/test_interface/test_lithium_ion.py
+++ b/tests/integration/test_models/test_submodels/test_interface/test_lithium_ion.py
@@ -97,7 +97,7 @@ def test_discretisation_lithium_ion(self):
 
         # Test
         whole_cell = ["negative electrode", "separator", "positive electrode"]
-        submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
         y = np.concatenate(
             [
                 submesh.nodes**2,
diff --git a/tests/integration/test_spatial_methods/test_finite_volume.py b/tests/integration/test_spatial_methods/test_finite_volume.py
index 67cf28ceaf..962e61804a 100644
--- a/tests/integration/test_spatial_methods/test_finite_volume.py
+++ b/tests/integration/test_spatial_methods/test_finite_volume.py
@@ -54,11 +54,11 @@ def get_error(n):
             disc.set_variable_slices([var])
 
             # Define exact solutions
-            combined_submesh = mesh.combine_submeshes(*whole_cell)
-            x = combined_submesh.nodes
+            submesh = mesh[whole_cell]
+            x = submesh.nodes
             y = np.sin(x) ** 2
             # var = sin(x)**2 --> dvardx = 2*sin(x)*cos(x)
-            x_edge = combined_submesh.edges
+            x_edge = submesh.edges
             grad_exact = 2 * np.sin(x_edge) * np.cos(x_edge)
 
             # Discretise and evaluate
@@ -90,8 +90,8 @@ def get_error(n):
             # create mesh and discretisation
             mesh = get_mesh_for_testing(n)
             disc = pybamm.Discretisation(mesh, spatial_methods)
-            combined_submesh = mesh.combine_submeshes(*whole_cell)
-            x = combined_submesh.nodes
+            submesh = mesh[whole_cell]
+            x = submesh.nodes
             x_edge = pybamm.standard_spatial_vars.x_edge
 
             # Define flux and eqn
diff --git a/tests/integration/test_spatial_methods/test_spectral_volume.py b/tests/integration/test_spatial_methods/test_spectral_volume.py
index a9e10ea78a..5a3ef6e88a 100644
--- a/tests/integration/test_spatial_methods/test_spectral_volume.py
+++ b/tests/integration/test_spatial_methods/test_spectral_volume.py
@@ -116,11 +116,11 @@ def get_error(n):
             disc.set_variable_slices([var])
 
             # Define exact solutions
-            combined_submesh = mesh.combine_submeshes(*whole_cell)
-            x = combined_submesh.nodes
+            submesh = mesh[whole_cell]
+            x = submesh.nodes
             y = np.sin(x) ** 2
             # var = sin(x)**2 --> dvardx = 2*sin(x)*cos(x)
-            x_edge = combined_submesh.edges
+            x_edge = submesh.edges
             grad_exact = 2 * np.sin(x_edge) * np.cos(x_edge)
 
             # Discretise and evaluate
@@ -154,8 +154,8 @@ def get_error(n):
             # create mesh and discretisation
             mesh = get_mesh_for_testing(n)
             disc = pybamm.Discretisation(mesh, spatial_methods)
-            combined_submesh = mesh.combine_submeshes(*whole_cell)
-            x = combined_submesh.nodes
+            submesh = mesh[whole_cell]
+            x = submesh.nodes
             x_edge = pybamm.standard_spatial_vars.x_edge
 
             # Define flux and bcs
diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 15f8c78646..14cd4f6503 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -258,9 +258,9 @@ def test_discretise_slicing(self):
 
         self.assertEqual(disc.y_slices, {c: [slice(0, 100)]})
 
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
-        c_true = combined_submesh.nodes**2
+        c_true = submesh.nodes**2
         y = c_true
         np.testing.assert_array_equal(y[disc.y_slices[c][0]], c_true)
 
@@ -280,7 +280,7 @@ def test_discretise_slicing(self):
         np.testing.assert_array_equal(
             disc.bounds[1], [np.inf] * 100 + [1] * 100 + [np.inf] * 40
         )
-        d_true = 4 * combined_submesh.nodes
+        d_true = 4 * submesh.nodes
         jn_true = mesh["negative electrode"].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
         np.testing.assert_array_equal(y[disc.y_slices[c][0]], c_true)
@@ -310,7 +310,7 @@ def test_discretise_slicing(self):
         np.testing.assert_array_equal(
             disc.bounds[1], [np.inf] * 100 + [1] * 100 + [np.inf] * 100
         )
-        d_true = 4 * combined_submesh.nodes
+        d_true = 4 * submesh.nodes
         jn_true = mesh["negative electrode"].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
         np.testing.assert_array_equal(y[disc.y_slices[c][0]], c_true)
@@ -453,8 +453,8 @@ def test_discretise_spatial_operator(self):
             self.assertIsInstance(eqn_disc, pybamm.MatrixMultiplication)
             self.assertIsInstance(eqn_disc.children[0], pybamm.Matrix)
 
-            combined_submesh = mesh.combine_submeshes(*whole_cell)
-            y = combined_submesh.nodes**2
+            submesh = mesh[whole_cell]
+            y = submesh.nodes**2
             var_disc = disc.process_symbol(var)
             # grad and var are identity operators here (for testing purposes)
             np.testing.assert_array_equal(
@@ -470,7 +470,7 @@ def test_discretise_spatial_operator(self):
             self.assertIsInstance(eqn_disc.children[1], pybamm.MatrixMultiplication)
             self.assertIsInstance(eqn_disc.children[1].children[0], pybamm.Matrix)
 
-            y = combined_submesh.nodes**2
+            y = submesh.nodes**2
             var_disc = disc.process_symbol(var)
             # grad and var are identity operators here (for testing purposes)
             np.testing.assert_array_equal(
@@ -511,9 +511,9 @@ def test_process_dict(self):
         disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
-        y = combined_submesh.nodes[:, np.newaxis] ** 2
+        y = submesh.nodes[:, np.newaxis] ** 2
         disc.bcs = boundary_conditions
 
         disc.set_variable_slices(list(rhs.keys()))
@@ -524,7 +524,7 @@ def test_process_dict(self):
         y0 = disc.process_dict(initial_conditions)
         np.testing.assert_array_equal(
             y0[c].evaluate(0, None),
-            3 * np.ones_like(combined_submesh.nodes[:, np.newaxis]),
+            3 * np.ones_like(submesh.nodes[:, np.newaxis]),
         )
         # vars
         processed_vars = disc.process_dict(variables)
@@ -538,9 +538,9 @@ def test_process_dict(self):
         rhs = {c: pybamm.div(N), T: pybamm.div(q)}
         initial_conditions = {c: pybamm.Scalar(3), T: pybamm.Scalar(5)}
         boundary_conditions = {}
-        y = np.concatenate(
-            [combined_submesh.nodes**2, mesh["negative electrode"].nodes ** 4]
-        )[:, np.newaxis]
+        y = np.concatenate([submesh.nodes**2, mesh["negative electrode"].nodes ** 4])[
+            :, np.newaxis
+        ]
 
         variables = list(rhs.keys())
         disc.set_variable_slices(variables)
@@ -556,7 +556,7 @@ def test_process_dict(self):
         y0 = disc.process_dict(initial_conditions)
         np.testing.assert_array_equal(
             y0[c].evaluate(0, None),
-            3 * np.ones_like(combined_submesh.nodes[:, np.newaxis]),
+            3 * np.ones_like(submesh.nodes[:, np.newaxis]),
         )
         np.testing.assert_array_equal(
             y0[T].evaluate(0, None),
@@ -589,12 +589,12 @@ def test_process_model_ode(self):
         disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
         disc.process_model(model)
 
         y0 = model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
-            y0, 3 * np.ones_like(combined_submesh.nodes[:, np.newaxis])
+            y0, 3 * np.ones_like(submesh.nodes[:, np.newaxis])
         )
         np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
 
@@ -604,15 +604,15 @@ def test_process_model_ode(self):
 
         # mass matrix is identity
         np.testing.assert_array_equal(
-            np.eye(combined_submesh.nodes.shape[0]), model.mass_matrix.entries.toarray()
+            np.eye(submesh.nodes.shape[0]), model.mass_matrix.entries.toarray()
         )
 
         # Create StateVector to differentiate model with respect to
-        y = pybamm.StateVector(slice(0, combined_submesh.npts))
+        y = pybamm.StateVector(slice(0, submesh.npts))
 
         # jacobian is identity
         jacobian = model.concatenated_rhs.jac(y).evaluate(0, y0)
-        np.testing.assert_array_equal(np.eye(combined_submesh.npts), jacobian.toarray())
+        np.testing.assert_array_equal(np.eye(submesh.npts), jacobian.toarray())
 
         # several equations
         T = pybamm.Variable("T", domain=["negative electrode"])
@@ -638,7 +638,7 @@ def test_process_model_ode(self):
         y0_expect = np.empty((0, 1))
         for var_id, _ in sorted(disc.y_slices.items(), key=lambda kv: kv[1]):
             if var_id == c:
-                vect = 2 * np.ones_like(combined_submesh.nodes[:, np.newaxis])
+                vect = 2 * np.ones_like(submesh.nodes[:, np.newaxis])
             elif var_id == T:
                 vect = 5 * np.ones_like(mesh["negative electrode"].nodes[:, np.newaxis])
             else:
@@ -744,36 +744,34 @@ def test_process_model_dae(self):
         mesh = disc.mesh
 
         disc.process_model(model)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         y0 = model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
             y0,
             np.concatenate(
                 [
-                    3 * np.ones_like(combined_submesh.nodes),
-                    6 * np.ones_like(combined_submesh.nodes),
+                    3 * np.ones_like(submesh.nodes),
+                    6 * np.ones_like(submesh.nodes),
                 ]
             )[:, np.newaxis],
         )
 
         # grad and div are identity operators here
         np.testing.assert_array_equal(
-            y0[: combined_submesh.npts], model.concatenated_rhs.evaluate(None, y0)
+            y0[: submesh.npts], model.concatenated_rhs.evaluate(None, y0)
         )
 
         np.testing.assert_array_equal(
             model.concatenated_algebraic.evaluate(None, y0),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
         # mass matrix is identity upper left, zeros elsewhere
         mass = block_diag(
             (
-                np.eye(np.size(combined_submesh.nodes)),
-                np.zeros(
-                    (np.size(combined_submesh.nodes), np.size(combined_submesh.nodes))
-                ),
+                np.eye(np.size(submesh.nodes)),
+                np.zeros((np.size(submesh.nodes), np.size(submesh.nodes))),
             )
         )
         np.testing.assert_array_equal(
@@ -789,17 +787,17 @@ def test_process_model_dae(self):
         jacobian_actual = np.block(
             [
                 [
-                    np.eye(np.size(combined_submesh.nodes)),
+                    np.eye(np.size(submesh.nodes)),
                     np.zeros(
                         (
-                            np.size(combined_submesh.nodes),
-                            np.size(combined_submesh.nodes),
+                            np.size(submesh.nodes),
+                            np.size(submesh.nodes),
                         )
                     ),
                 ],
                 [
-                    -2 * np.eye(np.size(combined_submesh.nodes)),
-                    np.eye(np.size(combined_submesh.nodes)),
+                    -2 * np.eye(np.size(submesh.nodes)),
+                    np.eye(np.size(submesh.nodes)),
                 ],
             ]
         )
@@ -840,12 +838,12 @@ def test_process_model_algebraic(self):
         mesh = disc.mesh
 
         disc.process_model(model)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         y0 = model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
             y0,
-            np.zeros_like(combined_submesh.nodes)[:, np.newaxis],
+            np.zeros_like(submesh.nodes)[:, np.newaxis],
         )
 
         # grad and div are identity operators here
@@ -855,19 +853,17 @@ def test_process_model_algebraic(self):
 
         np.testing.assert_array_equal(
             model.concatenated_algebraic.evaluate(None, y0),
-            -np.ones_like(combined_submesh.nodes[:, np.newaxis]),
+            -np.ones_like(submesh.nodes[:, np.newaxis]),
         )
 
         # mass matrix is identity upper left, zeros elsewhere
-        mass = np.zeros(
-            (np.size(combined_submesh.nodes), np.size(combined_submesh.nodes))
-        )
+        mass = np.zeros((np.size(submesh.nodes), np.size(submesh.nodes)))
         np.testing.assert_array_equal(mass, model.mass_matrix.entries.toarray())
 
         # jacobian
         y = pybamm.StateVector(slice(0, np.size(y0)))
         jacobian = model.concatenated_algebraic.jac(y).evaluate(0, y0)
-        np.testing.assert_array_equal(np.eye(combined_submesh.npts), jacobian.toarray())
+        np.testing.assert_array_equal(np.eye(submesh.npts), jacobian.toarray())
 
     def test_process_model_concatenation(self):
         # concatenation of variables as the key
@@ -889,14 +885,14 @@ def test_process_model_concatenation(self):
         disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        combined_submesh = mesh.combine_submeshes(
+        submesh = mesh[(
             "negative electrode", "separator", "positive electrode"
-        )
+        )]
 
         disc.process_model(model)
         y0 = model.concatenated_initial_conditions.evaluate()
         np.testing.assert_array_equal(
-            y0, 3 * np.ones_like(combined_submesh.nodes[:, np.newaxis])
+            y0, 3 * np.ones_like(submesh.nodes[:, np.newaxis])
         )
 
         # grad and div are identity operators here
@@ -961,13 +957,13 @@ def test_broadcast(self):
         disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # scalar
         broad = disc.process_symbol(pybamm.FullBroadcast(a, whole_cell, {}))
         np.testing.assert_array_equal(
             broad.evaluate(inputs={"a": 7}),
-            7 * np.ones_like(combined_submesh.nodes[:, np.newaxis]),
+            7 * np.ones_like(submesh.nodes[:, np.newaxis]),
         )
         self.assertEqual(broad.domain, whole_cell)
 
@@ -1312,7 +1308,7 @@ def test_process_input_variable(self):
 
         a = pybamm.InputParameter("a", ["negative electrode", "separator"])
         a_disc = disc.process_symbol(a)
-        n = disc.mesh.combine_submeshes(*a.domain).npts
+        n = disc.mesh[a.domain].npts
         self.assertEqual(a_disc._expected_size, n)
 
     def test_process_not_constant(self):
diff --git a/tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py b/tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py
index eee960e069..77bfecb51f 100644
--- a/tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py
+++ b/tests/unit/test_expression_tree/test_operations/test_evaluate_julia.py
@@ -233,8 +233,8 @@ def test_evaluator_julia_domain_concatenation(self):
         }
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
-        combined_submesh = mesh.combine_submeshes(*c.domain)
-        nodes = combined_submesh.nodes
+        submesh = mesh[c.domain]
+        nodes = submesh.nodes
         y_tests = [nodes**2 + 1, np.cos(nodes)]
 
         # discretise and evaluate the variable
@@ -264,10 +264,8 @@ def test_evaluator_julia_domain_concatenation_2D(self):
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
-        combined_submesh = mesh.combine_submeshes(*c.domain)
-        nodes = np.linspace(
-            0, 1, combined_submesh.npts * mesh["current collector"].npts
-        )
+        submesh = mesh[c.domain]
+        nodes = np.linspace(0, 1, submesh.npts * mesh["current collector"].npts)
         y_tests = [nodes**2 + 1, np.cos(nodes)]
 
         # discretise and evaluate the variable
@@ -283,7 +281,7 @@ def test_evaluator_julia_discretised_operators(self):
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         var = pybamm.Variable("var", domain=whole_cell)
         boundary_conditions = {
@@ -304,7 +302,7 @@ def test_evaluator_julia_discretised_operators(self):
         div_eqn_disc = disc.process_symbol(div_eqn)
 
         # test
-        nodes = combined_submesh.nodes
+        nodes = submesh.nodes
         y_tests = [nodes**2 + 1, np.cos(nodes)]
 
         for i, expr in enumerate([grad_eqn_disc, div_eqn_disc]):
diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py
index c42497f8f2..84aefebfec 100644
--- a/tests/unit/test_meshes/test_meshes.py
+++ b/tests/unit/test_meshes/test_meshes.py
@@ -82,6 +82,7 @@ def test_mesh_creation(self):
             mesh["positive electrode"].edges[0], mesh["separator"].edges[-1]
         )
         for domain in mesh:
+            domain = domain[0]
             if domain != "current collector":
                 self.assertEqual(len(mesh[domain].edges), len(mesh[domain].nodes) + 1)
 
@@ -223,7 +224,7 @@ def test_combine_submeshes(self):
         mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
 
         # create submesh
-        submesh = mesh.combine_submeshes("negative electrode", "separator")
+        submesh = mesh[("negative electrode", "separator")]
         self.assertEqual(submesh.edges[0], 0)
         self.assertEqual(submesh.edges[-1], mesh["separator"].edges[-1])
         np.testing.assert_almost_equal(
diff --git a/tests/unit/test_parameters/test_lead_acid_parameters.py b/tests/unit/test_parameters/test_lead_acid_parameters.py
index 49e9d5d716..bb67d74e21 100644
--- a/tests/unit/test_parameters/test_lead_acid_parameters.py
+++ b/tests/unit/test_parameters/test_lead_acid_parameters.py
@@ -67,10 +67,8 @@ def test_concatenated_parameters(self):
         processed_s = disc.process_symbol(parameter_values.process_symbol(s_param))
 
         # test output
-        combined_submeshes = disc.mesh.combine_submeshes(
-            "negative electrode", "separator", "positive electrode"
-        )
-        self.assertEqual(processed_s.shape, (combined_submeshes.npts, 1))
+        submeshes = disc.mesh[("negative electrode", "separator", "positive electrode")]
+        self.assertEqual(processed_s.shape, (submeshes.npts, 1))
 
     def test_current_functions(self):
         # create current functions
diff --git a/tests/unit/test_plotting/test_quick_plot.py b/tests/unit/test_plotting/test_quick_plot.py
index ff28a3fc48..a7816a748a 100644
--- a/tests/unit/test_plotting/test_quick_plot.py
+++ b/tests/unit/test_plotting/test_quick_plot.py
@@ -322,7 +322,7 @@ def test_loqs_spme(self):
             c_e_var = solution["Electrolyte concentration [mol.m-3]"]
             # 1D variables should be evaluated on edges
             L_x = param.evaluate(model.param.L_x)
-            c_e = c_e_var(t=t, x=mesh.combine_submeshes(*c_e_var.domain).edges * L_x)
+            c_e = c_e_var(t=t, x=mesh[c_e_var.domain].edges * L_x)
 
             for unit, scale in zip(["seconds", "minutes", "hours"], [1, 60, 3600]):
                 quick_plot = pybamm.QuickPlot(
diff --git a/tests/unit/test_solvers/test_scikits_solvers.py b/tests/unit/test_solvers/test_scikits_solvers.py
index 5ad45c2886..60aa60adf6 100644
--- a/tests/unit/test_solvers/test_scikits_solvers.py
+++ b/tests/unit/test_solvers/test_scikits_solvers.py
@@ -165,10 +165,8 @@ def test_model_solver_ode_jacobian_python(self):
 
         # Add user-supplied Jacobian to model
         mesh = get_mesh_for_testing()
-        combined_submesh = mesh.combine_submeshes(
-            "negative electrode", "separator", "positive electrode"
-        )
-        N = combined_submesh.npts
+        submesh = mesh[("negative electrode", "separator", "positive electrode")]
+        N = submesh.npts
 
         # Solve testing various linear solvers
         linsolvers = [
@@ -478,10 +476,8 @@ def test_model_solver_dae_with_jacobian_python(self):
 
         # Add user-supplied Jacobian to model
         mesh = get_mesh_for_testing()
-        combined_submesh = mesh.combine_submeshes(
-            "negative electrode", "separator", "positive electrode"
-        )
-        N = combined_submesh.npts
+        submesh = mesh[("negative electrode", "separator", "positive electrode")]
+        N = submesh.npts
 
         def jacobian(t, y):
             return np.block(
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index 63cd50a40d..c1de118ea3 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -123,10 +123,10 @@ def test_model_solver_ode_with_jacobian_python(self):
         disc.process_model(model)
 
         # Add user-supplied Jacobian to model
-        combined_submesh = mesh.combine_submeshes(
+        submesh = mesh[(
             "negative electrode", "separator", "positive electrode"
-        )
-        N = combined_submesh.npts
+        )]
+        N = submesh.npts
 
         # construct jacobian in order of model.rhs
         J = []
diff --git a/tests/unit/test_spatial_methods/test_base_spatial_method.py b/tests/unit/test_spatial_methods/test_base_spatial_method.py
index 259173092e..46a792f26a 100644
--- a/tests/unit/test_spatial_methods/test_base_spatial_method.py
+++ b/tests/unit/test_spatial_methods/test_base_spatial_method.py
@@ -107,7 +107,7 @@ def test_discretise_spatial_variable(self):
             var_disc = spatial_method.spatial_variable(var)
             self.assertIsInstance(var_disc, pybamm.Vector)
             np.testing.assert_array_equal(
-                var_disc.evaluate()[:, 0], mesh.combine_submeshes(*var.domain).nodes
+                var_disc.evaluate()[:, 0], mesh[var.domain].nodes
             )
 
         # edges
@@ -118,7 +118,7 @@ def test_discretise_spatial_variable(self):
             var_disc = spatial_method.spatial_variable(var)
             self.assertIsInstance(var_disc, pybamm.Vector)
             np.testing.assert_array_equal(
-                var_disc.evaluate()[:, 0], mesh.combine_submeshes(*var.domain).edges
+                var_disc.evaluate()[:, 0], mesh[var.domain].edges
             )
 
     def test_boundary_value_checks(self):
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py b/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py
index 07462cf3d2..2521755708 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_extrapolation.py
@@ -248,7 +248,7 @@ def test_linear_extrapolate_left_right(self):
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
-        macro_submesh = mesh.combine_submeshes(*whole_cell)
+        macro_submesh = mesh[whole_cell]
         micro_submesh = mesh["negative particle"]
 
         # Macroscale
@@ -315,7 +315,7 @@ def test_quadratic_extrapolate_left_right(self):
         disc = pybamm.Discretisation(mesh, spatial_methods)
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
-        macro_submesh = mesh.combine_submeshes(*whole_cell)
+        macro_submesh = mesh[whole_cell]
         micro_submesh = mesh["negative particle"]
 
         # Macroscale
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
index 008e39b26c..87a4321f9b 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_finite_volume.py
@@ -64,9 +64,7 @@ def test_concatenation(self):
         edges = [pybamm.Vector(mesh[dom].edges, domain=dom) for dom in whole_cell]
         # Concatenation of edges should get averaged to nodes first, using edge_to_node
         v_disc = fin_vol.concatenation(edges)
-        np.testing.assert_array_equal(
-            v_disc.evaluate()[:, 0], mesh.combine_submeshes(*whole_cell).nodes
-        )
+        np.testing.assert_array_equal(v_disc.evaluate()[:, 0], mesh[whole_cell].nodes)
 
         # test for bad shape
         edges = [
@@ -81,12 +79,12 @@ def test_discretise_diffusivity_times_spatial_operator(self):
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
         whole_cell = ["negative electrode", "separator", "positive electrode"]
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # Discretise some equations where averaging is needed
         var = pybamm.Variable("var", domain=whole_cell)
         disc.set_variable_slices([var])
-        y_test = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        y_test = np.ones_like(submesh.nodes[:, np.newaxis])
         for eqn in [
             var * pybamm.grad(var),
             var**2 * pybamm.grad(var),
@@ -165,9 +163,7 @@ def test_discretise_spatial_variable(self):
         self.assertIsInstance(x2_disc, pybamm.Vector)
         np.testing.assert_array_equal(
             x2_disc.evaluate(),
-            disc.mesh.combine_submeshes("negative electrode", "separator").nodes[
-                :, np.newaxis
-            ],
+            disc.mesh[("negative electrode", "separator")].nodes[:, np.newaxis],
         )
         # microscale
         r = 3 * pybamm.SpatialVariable("r", ["negative particle"])
@@ -194,11 +190,11 @@ def test_mass_matrix_shape(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
         disc.process_model(model)
 
         # Mass matrix
-        mass = np.eye(combined_submesh.npts)
+        mass = np.eye(submesh.npts)
         np.testing.assert_array_equal(mass, model.mass_matrix.entries.toarray())
 
     def test_p2d_mass_matrix_shape(self):
@@ -238,15 +234,15 @@ def test_jacobian(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
         spatial_method = pybamm.FiniteVolume()
         spatial_method.build(mesh)
 
         # Setup variable
         var = pybamm.Variable("var", domain=whole_cell)
         disc.set_variable_slices([var])
-        y = pybamm.StateVector(slice(0, combined_submesh.npts))
-        y_test = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        y = pybamm.StateVector(slice(0, submesh.npts))
+        y_test = np.ones_like(submesh.nodes[:, np.newaxis])
 
         # grad
         eqn = pybamm.grad(var)
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py b/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py
index d9f35523f7..0a9d1dd102 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_ghost_nodes_and_neumann.py
@@ -30,8 +30,8 @@ def test_add_ghost_nodes(self):
         sp_meth = pybamm.FiniteVolume()
         sp_meth.build(mesh)
         sym_ghost, _ = sp_meth.add_ghost_nodes(var, discretised_symbol, bcs)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
-        y_test = np.linspace(0, 1, combined_submesh.npts)
+        submesh = mesh[whole_cell]
+        y_test = np.linspace(0, 1, submesh.npts)
         np.testing.assert_array_equal(
             sym_ghost.evaluate(y=y_test)[1:-1], discretised_symbol.evaluate(y=y_test)
         )
@@ -76,8 +76,8 @@ def test_add_ghost_nodes_concatenation(self):
         }
 
         # Test
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
-        y_test = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        submesh = mesh[whole_cell]
+        y_test = np.ones_like(submesh.nodes[:, np.newaxis])
 
         # both
         sp_meth = pybamm.FiniteVolume()
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py b/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
index a586939a08..ac5ec6f489 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_grad_div_shapes.py
@@ -22,11 +22,11 @@ def test_grad_div_shapes_Dirichlet_bcs(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # Test gradient of constant is zero
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
         boundary_conditions = {
@@ -40,11 +40,11 @@ def test_grad_div_shapes_Dirichlet_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test operations on linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -58,13 +58,13 @@ def test_grad_div_shapes_Dirichlet_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_cylindrical_grad_div_shapes_Dirichlet_bcs(self):
@@ -261,11 +261,11 @@ def test_grad_div_shapes_Neumann_bcs(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # Test gradient of constant is zero
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
         boundary_conditions = {
@@ -279,11 +279,11 @@ def test_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test operations on linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -297,13 +297,13 @@ def test_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
@@ -316,10 +316,10 @@ def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # Test gradient and divergence of a constant
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
         N = pybamm.grad(var)
@@ -336,17 +336,17 @@ def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(1)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
         # Test gradient and divergence of linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(1), "Neumann"),
@@ -358,13 +358,13 @@ def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_cylindrical_grad_div_shapes_Neumann_bcs(self):
@@ -439,13 +439,13 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"negative particle": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes("negative particle")
+        submesh = mesh["negative particle"]
 
         # Test gradient
         var = pybamm.Variable("var", domain="negative particle")
         grad_eqn = pybamm.grad(var)
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(0), "Neumann"),
@@ -457,10 +457,10 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
         # grad(r) == 1
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(1), "Neumann"),
@@ -471,12 +471,12 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test divergence of gradient
         # div(grad(r^2)) = 6 , N_left = 0, N_right = 2
-        quadratic_y = combined_submesh.nodes**2
+        quadratic_y = submesh.nodes**2
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -489,7 +489,7 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, quadratic_y),
-            6 * np.ones((combined_submesh.npts, 1)),
+            6 * np.ones((submesh.npts, 1)),
         )
 
     def test_p2d_spherical_grad_div_shapes_Neumann_bcs(self):
@@ -548,11 +548,11 @@ def test_grad_div_shapes_mixed_domain(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.FiniteVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        submesh = mesh[("negative electrode", "separator")]
 
         # Test gradient of constant
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=["negative electrode", "separator"])
         grad_eqn = pybamm.grad(var)
         boundary_conditions = {
@@ -566,17 +566,17 @@ def test_grad_div_shapes_mixed_domain(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test operations on linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(0), "Dirichlet"),
-                "right": (pybamm.Scalar(combined_submesh.edges[-1]), "Dirichlet"),
+                "right": (pybamm.Scalar(submesh.edges[-1]), "Dirichlet"),
             }
         }
         disc.bcs = boundary_conditions
@@ -584,13 +584,13 @@ def test_grad_div_shapes_mixed_domain(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_grad_1plus1d(self):
@@ -627,13 +627,11 @@ def test_grad_1plus1d(self):
         grad_eqn_disc = disc.process_symbol(pybamm.grad(var))
 
         # Evaulate
-        combined_submesh = mesh.combine_submeshes(*var.domain)
-        linear_y = np.outer(np.linspace(0, 1, 15), combined_submesh.nodes).reshape(
+        submesh = mesh[var.domain]
+        linear_y = np.outer(np.linspace(0, 1, 15), submesh.nodes).reshape(-1, 1)
+        expected = np.outer(np.linspace(0, 1, 15), np.ones_like(submesh.edges)).reshape(
             -1, 1
         )
-        expected = np.outer(
-            np.linspace(0, 1, 15), np.ones_like(combined_submesh.edges)
-        ).reshape(-1, 1)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y), expected
         )
diff --git a/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py b/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py
index 0885eb58ca..8b021751c1 100644
--- a/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py
+++ b/tests/unit/test_spatial_methods/test_finite_volume/test_integration.py
@@ -33,15 +33,15 @@ def test_definite_integral(self):
         integral_eqn = pybamm.Integral(var, x)
         disc.set_variable_slices([var])
         integral_eqn_disc = disc.process_symbol(integral_eqn)
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        submesh = mesh[("negative electrode", "separator")]
 
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), ln + ls)
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         np.testing.assert_array_almost_equal(
             integral_eqn_disc.evaluate(None, linear_y), (ln + ls) ** 2 / 2
         )
-        cos_y = np.cos(combined_submesh.nodes[:, np.newaxis])
+        cos_y = np.cos(submesh.nodes[:, np.newaxis])
         np.testing.assert_array_almost_equal(
             integral_eqn_disc.evaluate(None, cos_y), np.sin(ln + ls), decimal=4
         )
@@ -52,15 +52,15 @@ def test_definite_integral(self):
         integral_eqn = pybamm.Integral(var, x)
         disc.set_variable_slices([var])
         integral_eqn_disc = disc.process_symbol(integral_eqn)
-        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+        submesh = mesh[("separator", "positive electrode")]
 
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), ls + lp)
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         self.assertAlmostEqual(
             integral_eqn_disc.evaluate(None, linear_y)[0][0], (1 - (ln) ** 2) / 2
         )
-        cos_y = np.cos(combined_submesh.nodes[:, np.newaxis])
+        cos_y = np.cos(submesh.nodes[:, np.newaxis])
         np.testing.assert_array_almost_equal(
             integral_eqn_disc.evaluate(None, cos_y), np.sin(1) - np.sin(ln), decimal=4
         )
@@ -328,22 +328,22 @@ def test_indefinite_integral(self):
         left_boundary_value = pybamm.BoundaryValue(int_grad_phi, "left")
         left_boundary_value_disc = disc.process_symbol(left_boundary_value)
 
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        submesh = mesh[("negative electrode", "separator")]
 
         # constant case
-        phi_exact = np.ones((combined_submesh.npts, 1))
+        phi_exact = np.ones((submesh.npts, 1))
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         phi_approx += 1  # add constant of integration
         np.testing.assert_array_equal(phi_exact, phi_approx)
         self.assertEqual(left_boundary_value_disc.evaluate(y=phi_exact), 0)
         # linear case
-        phi_exact = combined_submesh.nodes[:, np.newaxis]
+        phi_exact = submesh.nodes[:, np.newaxis]
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(phi_exact, phi_approx)
         self.assertEqual(left_boundary_value_disc.evaluate(y=phi_exact), 0)
 
         # sine case
-        phi_exact = np.sin(combined_submesh.nodes[:, np.newaxis])
+        phi_exact = np.sin(submesh.nodes[:, np.newaxis])
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(phi_exact, phi_approx)
         self.assertEqual(left_boundary_value_disc.evaluate(y=phi_exact), 0)
@@ -364,17 +364,17 @@ def test_indefinite_integral(self):
         int_grad_phi_disc = disc.process_symbol(int_grad_phi)
         left_boundary_value = pybamm.BoundaryValue(int_grad_phi, "left")
         left_boundary_value_disc = disc.process_symbol(left_boundary_value)
-        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+        submesh = mesh[("separator", "positive electrode")]
 
         # constant case
-        phi_exact = np.ones((combined_submesh.npts, 1))
+        phi_exact = np.ones((submesh.npts, 1))
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         phi_approx += 1  # add constant of integration
         np.testing.assert_array_equal(phi_exact, phi_approx)
         self.assertEqual(left_boundary_value_disc.evaluate(y=phi_exact), 0)
 
         # linear case
-        phi_exact = combined_submesh.nodes[:, np.newaxis] - combined_submesh.edges[0]
+        phi_exact = submesh.nodes[:, np.newaxis] - submesh.edges[0]
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(phi_exact, phi_approx)
         np.testing.assert_array_almost_equal(
@@ -382,9 +382,7 @@ def test_indefinite_integral(self):
         )
 
         # sine case
-        phi_exact = np.sin(
-            combined_submesh.nodes[:, np.newaxis] - combined_submesh.edges[0]
-        )
+        phi_exact = np.sin(submesh.nodes[:, np.newaxis] - submesh.edges[0])
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(phi_exact, phi_approx)
         np.testing.assert_array_almost_equal(
@@ -427,17 +425,17 @@ def test_indefinite_integral(self):
         c_integral_disc = disc.process_symbol(c_integral)
         left_boundary_value = pybamm.BoundaryValue(c_integral, "left")
         left_boundary_value_disc = disc.process_symbol(left_boundary_value)
-        combined_submesh = mesh["negative particle"]
+        submesh = mesh["negative particle"]
 
         # constant case
-        c_exact = np.ones((combined_submesh.npts, 1))
+        c_exact = np.ones((submesh.npts, 1))
         c_approx = c_integral_disc.evaluate(None, c_exact)
         c_approx += 1  # add constant of integration
         np.testing.assert_array_equal(c_exact, c_approx)
         self.assertEqual(left_boundary_value_disc.evaluate(y=c_exact), 0)
 
         # linear case
-        c_exact = combined_submesh.nodes[:, np.newaxis]
+        c_exact = submesh.nodes[:, np.newaxis]
         c_approx = c_integral_disc.evaluate(None, c_exact)
         np.testing.assert_array_almost_equal(c_exact, c_approx)
         np.testing.assert_array_almost_equal(
@@ -445,7 +443,7 @@ def test_indefinite_integral(self):
         )
 
         # sine case
-        c_exact = np.sin(combined_submesh.nodes[:, np.newaxis])
+        c_exact = np.sin(submesh.nodes[:, np.newaxis])
         c_approx = c_integral_disc.evaluate(None, c_exact)
         np.testing.assert_array_almost_equal(c_exact, c_approx, decimal=3)
         np.testing.assert_array_almost_equal(
@@ -475,18 +473,18 @@ def test_backward_indefinite_integral(self):
         int_grad_phi_disc = disc.process_symbol(int_grad_phi)
         right_boundary_value = pybamm.BoundaryValue(int_grad_phi, "right")
         right_boundary_value_disc = disc.process_symbol(right_boundary_value)
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        submesh = mesh[("negative electrode", "separator")]
 
         # Test that the backward_integral(grad(phi)) = -phi
         # constant case
-        phi_exact = np.ones((combined_submesh.npts, 1))
+        phi_exact = np.ones((submesh.npts, 1))
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact)
         phi_approx += 1  # add constant of integration
         np.testing.assert_array_equal(phi_exact, phi_approx)
         self.assertEqual(right_boundary_value_disc.evaluate(y=phi_exact), 0)
 
         # linear case
-        phi_exact = combined_submesh.nodes - combined_submesh.edges[-1]
+        phi_exact = submesh.nodes - submesh.edges[-1]
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(phi_exact, -phi_approx)
         np.testing.assert_array_almost_equal(
@@ -494,7 +492,7 @@ def test_backward_indefinite_integral(self):
         )
 
         # sine case
-        phi_exact = np.sin(combined_submesh.nodes - combined_submesh.edges[-1])
+        phi_exact = np.sin(submesh.nodes - submesh.edges[-1])
         phi_approx = int_grad_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(phi_exact, -phi_approx)
         np.testing.assert_array_almost_equal(
@@ -518,8 +516,8 @@ def test_indefinite_integral_of_broadcasted_to_cell_edges(self):
         i = pybamm.PrimaryBroadcastToEdges(1, phi.domain)
         x = pybamm.SpatialVariable("x", phi.domain)
         disc.set_variable_slices([phi])
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
-        x_end = combined_submesh.edges[-1]
+        submesh = mesh[("negative electrode", "separator")]
+        x_end = submesh.edges[-1]
 
         # take indefinite integral
         int_phi = pybamm.IndefiniteIntegral(i * phi, x)
@@ -528,12 +526,12 @@ def test_indefinite_integral_of_broadcasted_to_cell_edges(self):
         int_int_phi_disc = disc.process_symbol(int_int_phi)
 
         # constant case
-        phi_exact = np.ones_like(combined_submesh.nodes)
+        phi_exact = np.ones_like(submesh.nodes)
         phi_approx = int_int_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(x_end**2 / 2, phi_approx)
 
         # linear case
-        phi_exact = combined_submesh.nodes[:, np.newaxis]
+        phi_exact = submesh.nodes[:, np.newaxis]
         phi_approx = int_int_phi_disc.evaluate(None, phi_exact)
         np.testing.assert_array_almost_equal(x_end**3 / 6, phi_approx, decimal=4)
 
@@ -549,21 +547,21 @@ def test_indefinite_integral_on_nodes(self):
         disc.set_variable_slices([phi])
         int_phi_disc = disc.process_symbol(int_phi)
 
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        submesh = mesh[("negative electrode", "separator")]
 
         # constant case
-        phi_exact = np.ones((combined_submesh.npts, 1))
-        int_phi_exact = combined_submesh.edges
+        phi_exact = np.ones((submesh.npts, 1))
+        int_phi_exact = submesh.edges
         int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_equal(int_phi_exact, int_phi_approx)
         # linear case
-        phi_exact = combined_submesh.nodes
-        int_phi_exact = combined_submesh.edges**2 / 2
+        phi_exact = submesh.nodes
+        int_phi_exact = submesh.edges**2 / 2
         int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(int_phi_exact, int_phi_approx)
         # cos case
-        phi_exact = np.cos(combined_submesh.nodes)
-        int_phi_exact = np.sin(combined_submesh.edges)
+        phi_exact = np.cos(submesh.nodes)
+        int_phi_exact = np.sin(submesh.edges)
         int_phi_approx = int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(int_phi_exact, int_phi_approx, decimal=5)
 
@@ -595,21 +593,21 @@ def test_backward_indefinite_integral_on_nodes(self):
         disc.set_variable_slices([phi])
         back_int_phi_disc = disc.process_symbol(back_int_phi)
 
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
-        edges = combined_submesh.edges
+        submesh = mesh[("negative electrode", "separator")]
+        edges = submesh.edges
 
         # constant case
-        phi_exact = np.ones((combined_submesh.npts, 1))
+        phi_exact = np.ones((submesh.npts, 1))
         back_int_phi_exact = edges[-1] - edges
         back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(back_int_phi_exact, back_int_phi_approx)
         # linear case
-        phi_exact = combined_submesh.nodes
+        phi_exact = submesh.nodes
         back_int_phi_exact = edges[-1] ** 2 / 2 - edges**2 / 2
         back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(back_int_phi_exact, back_int_phi_approx)
         # cos case
-        phi_exact = np.cos(combined_submesh.nodes)
+        phi_exact = np.cos(submesh.nodes)
         back_int_phi_exact = np.sin(edges[-1]) - np.sin(edges)
         back_int_phi_approx = back_int_phi_disc.evaluate(None, phi_exact).flatten()
         np.testing.assert_array_almost_equal(
@@ -637,13 +635,13 @@ def test_forward_plus_backward_integral(self):
         ) + pybamm.BackwardIndefiniteIntegral(phi, x)
         int_plus_back_int_phi_disc = disc.process_symbol(int_plus_back_int_phi)
 
-        combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
+        submesh = mesh[("separator", "positive electrode")]
 
         # test
         for phi_exact in [
-            np.ones((combined_submesh.npts, 1)),
-            combined_submesh.nodes,
-            np.cos(combined_submesh.nodes),
+            np.ones((submesh.npts, 1)),
+            submesh.nodes,
+            np.cos(submesh.nodes),
         ]:
             np.testing.assert_array_almost_equal(
                 full_int_phi_disc.evaluate(y=phi_exact).flatten(),
diff --git a/tests/unit/test_spatial_methods/test_spectral_volume.py b/tests/unit/test_spatial_methods/test_spectral_volume.py
index ee63fc16c7..c2e5f3396d 100644
--- a/tests/unit/test_spatial_methods/test_spectral_volume.py
+++ b/tests/unit/test_spatial_methods/test_spectral_volume.py
@@ -121,11 +121,11 @@ def test_grad_div_shapes_Dirichlet_bcs(self):
         mesh = get_mesh_for_testing(1)
         spatial_methods = {"macroscale": pybamm.SpectralVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # Test gradient of constant is zero
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
         boundary_conditions = {
@@ -139,11 +139,11 @@ def test_grad_div_shapes_Dirichlet_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test operations on linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -157,13 +157,13 @@ def test_grad_div_shapes_Dirichlet_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_spherical_grad_div_shapes_Dirichlet_bcs(self):
@@ -312,11 +312,11 @@ def test_grad_div_shapes_Neumann_bcs(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.SpectralVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # Test gradient of constant is zero
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
         boundary_conditions = {
@@ -330,11 +330,11 @@ def test_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test operations on linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -348,13 +348,13 @@ def test_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
@@ -367,10 +367,10 @@ def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.SpectralVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes(*whole_cell)
+        submesh = mesh[whole_cell]
 
         # Test gradient and divergence of a constant
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
         N = pybamm.grad(var)
@@ -387,17 +387,17 @@ def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(1)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
         # Test gradient and divergence of linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(1), "Neumann"),
@@ -409,13 +409,13 @@ def test_grad_div_shapes_Dirichlet_and_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_spherical_grad_div_shapes_Neumann_bcs(self):
@@ -427,13 +427,13 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"negative particle": pybamm.SpectralVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes("negative particle")
+        submesh = mesh["negative particle"]
 
         # Test gradient
         var = pybamm.Variable("var", domain="negative particle")
         grad_eqn = pybamm.grad(var)
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(0), "Neumann"),
@@ -445,10 +445,10 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
         # grad(r) == 1
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(1), "Neumann"),
@@ -459,12 +459,12 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test divergence of gradient
         # div(grad(r^2)) = 6 , N_left = 0, N_right = 2
-        quadratic_y = combined_submesh.nodes**2
+        quadratic_y = submesh.nodes**2
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -477,7 +477,7 @@ def test_spherical_grad_div_shapes_Neumann_bcs(self):
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, quadratic_y),
-            6 * np.ones((combined_submesh.npts, 1)),
+            6 * np.ones((submesh.npts, 1)),
         )
 
     def test_p2d_spherical_grad_div_shapes_Neumann_bcs(self):
@@ -537,11 +537,11 @@ def test_grad_div_shapes_mixed_domain(self):
         mesh = get_mesh_for_testing()
         spatial_methods = {"macroscale": pybamm.SpectralVolume()}
         disc = pybamm.Discretisation(mesh, spatial_methods)
-        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        submesh = mesh[("negative electrode", "separator")]
 
         # Test gradient of constant
         # grad(1) = 0
-        constant_y = np.ones_like(combined_submesh.nodes[:, np.newaxis])
+        constant_y = np.ones_like(submesh.nodes[:, np.newaxis])
         var = pybamm.Variable("var", domain=["negative electrode", "separator"])
         grad_eqn = pybamm.grad(var)
         boundary_conditions = {
@@ -555,17 +555,17 @@ def test_grad_div_shapes_mixed_domain(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[:, np.newaxis]),
+            np.zeros_like(submesh.edges[:, np.newaxis]),
         )
 
         # Test operations on linear x
-        linear_y = combined_submesh.nodes
+        linear_y = submesh.nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
             var: {
                 "left": (pybamm.Scalar(0), "Dirichlet"),
-                "right": (pybamm.Scalar(combined_submesh.edges[-1]), "Dirichlet"),
+                "right": (pybamm.Scalar(submesh.edges[-1]), "Dirichlet"),
             }
         }
         disc.bcs = boundary_conditions
@@ -573,13 +573,13 @@ def test_grad_div_shapes_mixed_domain(self):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[:, np.newaxis]),
+            np.ones_like(submesh.edges[:, np.newaxis]),
         )
         # div(grad(x)) = 0
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
             div_eqn_disc.evaluate(None, linear_y),
-            np.zeros_like(combined_submesh.nodes[:, np.newaxis]),
+            np.zeros_like(submesh.nodes[:, np.newaxis]),
         )
 
     def test_grad_1plus1d(self):
@@ -616,13 +616,11 @@ def test_grad_1plus1d(self):
         grad_eqn_disc = disc.process_symbol(pybamm.grad(var))
 
         # Evaulate
-        combined_submesh = mesh.combine_submeshes(*var.domain)
-        linear_y = np.outer(np.linspace(0, 1, 15), combined_submesh.nodes).reshape(
+        submesh = mesh[var.domain]
+        linear_y = np.outer(np.linspace(0, 1, 15), submesh.nodes).reshape(-1, 1)
+        expected = np.outer(np.linspace(0, 1, 15), np.ones_like(submesh.edges)).reshape(
             -1, 1
         )
-        expected = np.outer(
-            np.linspace(0, 1, 15), np.ones_like(combined_submesh.edges)
-        ).reshape(-1, 1)
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y), expected
         )
diff --git a/tests/unit/test_spatial_methods/test_zero_dimensional_method.py b/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
index 4ba2e71cc6..b853080791 100644
--- a/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
+++ b/tests/unit/test_spatial_methods/test_zero_dimensional_method.py
@@ -39,7 +39,7 @@ def test_discretise_spatial_variable(self):
             var_disc = spatial_method.spatial_variable(var)
             self.assertIsInstance(var_disc, pybamm.Vector)
             np.testing.assert_array_equal(
-                var_disc.evaluate()[:, 0], mesh.combine_submeshes(*var.domain).nodes
+                var_disc.evaluate()[:, 0], mesh[var.domain].nodes
             )
 
         # edges
@@ -50,7 +50,7 @@ def test_discretise_spatial_variable(self):
             var_disc = spatial_method.spatial_variable(var)
             self.assertIsInstance(var_disc, pybamm.Vector)
             np.testing.assert_array_equal(
-                var_disc.evaluate()[:, 0], mesh.combine_submeshes(*var.domain).edges
+                var_disc.evaluate()[:, 0], mesh[var.domain].edges
             )
 
     def test_averages(self):

From 5d8196aaa30b22d1a359f0a9f440518cfd01d72b Mon Sep 17 00:00:00 2001
From: "pre-commit-ci[bot]"
 <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Date: Wed, 9 Nov 2022 19:03:56 +0000
Subject: [PATCH 06/12] style: pre-commit fixes

---
 tests/unit/test_discretisations/test_discretisation.py | 4 +---
 tests/unit/test_solvers/test_scipy_solver.py           | 4 +---
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/tests/unit/test_discretisations/test_discretisation.py b/tests/unit/test_discretisations/test_discretisation.py
index 14cd4f6503..01dd02a2b4 100644
--- a/tests/unit/test_discretisations/test_discretisation.py
+++ b/tests/unit/test_discretisations/test_discretisation.py
@@ -885,9 +885,7 @@ def test_process_model_concatenation(self):
         disc = get_discretisation_for_testing()
         mesh = disc.mesh
 
-        submesh = mesh[(
-            "negative electrode", "separator", "positive electrode"
-        )]
+        submesh = mesh[("negative electrode", "separator", "positive electrode")]
 
         disc.process_model(model)
         y0 = model.concatenated_initial_conditions.evaluate()
diff --git a/tests/unit/test_solvers/test_scipy_solver.py b/tests/unit/test_solvers/test_scipy_solver.py
index c1de118ea3..5603bd935a 100644
--- a/tests/unit/test_solvers/test_scipy_solver.py
+++ b/tests/unit/test_solvers/test_scipy_solver.py
@@ -123,9 +123,7 @@ def test_model_solver_ode_with_jacobian_python(self):
         disc.process_model(model)
 
         # Add user-supplied Jacobian to model
-        submesh = mesh[(
-            "negative electrode", "separator", "positive electrode"
-        )]
+        submesh = mesh[("negative electrode", "separator", "positive electrode")]
         N = submesh.npts
 
         # construct jacobian in order of model.rhs

From 8511bcc70d0e48ee8fdd34e518341331bbe3e0ed Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 14:05:54 -0500
Subject: [PATCH 07/12] flake8

---
 examples/scripts/compare_comsol/compare_comsol_DFN.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/scripts/compare_comsol/compare_comsol_DFN.py b/examples/scripts/compare_comsol/compare_comsol_DFN.py
index e5a9930a56..1b0497eecd 100644
--- a/examples/scripts/compare_comsol/compare_comsol_DFN.py
+++ b/examples/scripts/compare_comsol/compare_comsol_DFN.py
@@ -88,7 +88,7 @@ def get_interp_fun(variable_name, domain):
     )
 
     fun.domains = {"primary": domain}
-    fun.mesh = mesh.[domain]
+    fun.mesh = mesh[domain]
     fun.secondary_mesh = None
     return fun
 

From 7370e7002dd9e3cd362d7ef1492d78acbbe41fce Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 14:08:08 -0500
Subject: [PATCH 08/12] changelog

---
 CHANGELOG.md                               | 1 +
 pybamm/expression_tree/binary_operators.py | 9 ---------
 2 files changed, 1 insertion(+), 9 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 29d047d748..3926682e27 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@
 
 ## Optimizations
 
+- Added more rules for simplifying expressions, especially around Concatenations. Also, meshes constructed from multiple domains are now memoized ([#2443](https://github.com/pybamm-team/PyBaMM/pull/2443))
 - Added more rules for simplifying expressions. Constants in binary operators are now moved to the left by default (e.g. `x*2` returns `2*x`) ([#2424](https://github.com/pybamm-team/PyBaMM/pull/2424))
 
 ## Breaking changes
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 5712887773..d2e557e401 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -1049,15 +1049,6 @@ def simplified_multiplication(left, right):
                     elif isinstance(right, Subtraction):
                         return (left * r_left) - (left * r_right)
 
-    # # Move constants in multiplications to the left
-    # if isinstance(left, Multiplication) and left.left.is_constant():
-    #     l_left, l_right = left.orphans
-    #     print(l_left, l_right, right)
-    #     return l_left * (l_right * right)
-    # elif isinstance(right, Multiplication) and right.left.is_constant():
-    #     r_left, r_right = right.orphans
-    #     return r_left * (left * r_right)
-
     # Cancelling out common terms
     if isinstance(left, Division):
         # Simplify (a / b) * b to a

From 002657fe25c90380bca4afcaf726ff119908fafd Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 14:37:44 -0500
Subject: [PATCH 09/12] fix test

---
 pybamm/meshes/meshes.py                           | 2 ++
 tests/unit/test_meshes/test_meshes.py             | 3 +--
 tests/unit/test_meshes/test_scikit_fem_submesh.py | 8 ++++----
 3 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 9f5208e2ab..5e1bb0cd7a 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -111,8 +111,10 @@ def __init__(self, geometry, submesh_types, var_pts):
                         geometry[domain][spatial_variable][lim] = sym_eval
 
         # Create submeshes
+        self.base_domains = []
         for domain in geometry:
             self[domain] = submesh_types[domain](geometry[domain], submesh_pts[domain])
+            self.base_domains.append(domain)
 
         # add ghost meshes
         self.add_ghost_meshes()
diff --git a/tests/unit/test_meshes/test_meshes.py b/tests/unit/test_meshes/test_meshes.py
index 84aefebfec..b04402128b 100644
--- a/tests/unit/test_meshes/test_meshes.py
+++ b/tests/unit/test_meshes/test_meshes.py
@@ -81,8 +81,7 @@ def test_mesh_creation(self):
         self.assertEqual(
             mesh["positive electrode"].edges[0], mesh["separator"].edges[-1]
         )
-        for domain in mesh:
-            domain = domain[0]
+        for domain in mesh.base_domains:
             if domain != "current collector":
                 self.assertEqual(len(mesh[domain].edges), len(mesh[domain].nodes) + 1)
 
diff --git a/tests/unit/test_meshes/test_scikit_fem_submesh.py b/tests/unit/test_meshes/test_scikit_fem_submesh.py
index 91c2899da8..0e9c5561d4 100644
--- a/tests/unit/test_meshes/test_scikit_fem_submesh.py
+++ b/tests/unit/test_meshes/test_scikit_fem_submesh.py
@@ -52,7 +52,7 @@ def test_mesh_creation(self):
         self.assertEqual(
             mesh["positive electrode"].edges[0], mesh["separator"].edges[-1]
         )
-        for domain in mesh:
+        for domain in mesh.base_domains:
             if domain == "current collector":
                 # NOTE: only for degree 1
                 npts = var_pts["y"] * var_pts["z"]
@@ -223,7 +223,7 @@ def test_mesh_creation(self):
         self.assertEqual(
             mesh["positive electrode"].edges[0], mesh["separator"].edges[-1]
         )
-        for domain in mesh:
+        for domain in mesh.base_domains:
             if domain == "current collector":
                 # NOTE: only for degree 1
                 npts = var_pts["y"] * var_pts["z"]
@@ -303,7 +303,7 @@ def test_mesh_creation(self):
         self.assertEqual(
             mesh["positive electrode"].edges[0], mesh["separator"].edges[-1]
         )
-        for domain in mesh:
+        for domain in mesh.base_domains:
             if domain == "current collector":
                 # NOTE: only for degree 1
                 npts = var_pts["y"] * var_pts["z"]
@@ -391,7 +391,7 @@ def test_mesh_creation(self):
         self.assertEqual(
             mesh["positive electrode"].edges[0], mesh["separator"].edges[-1]
         )
-        for domain in mesh:
+        for domain in mesh.base_domains:
             if domain == "current collector":
                 # NOTE: only for degree 1
                 npts = var_pts["y"] * var_pts["z"]

From 96e1cb25f2b0c29d18f51a2c6097c012c9678e42 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 15:11:01 -0500
Subject: [PATCH 10/12] fix integration test

---
 tests/integration/test_models/standard_output_tests.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py
index 295837a63c..7add30a0e5 100644
--- a/tests/integration/test_models/standard_output_tests.py
+++ b/tests/integration/test_models/standard_output_tests.py
@@ -550,7 +550,7 @@ def test_splitting(self):
             (self.c_e_n(t, x_n), self.c_e_s(t, x_s), self.c_e_p(t, x_p)), axis=0
         )
 
-        np.testing.assert_array_equal(self.c_e(t, x), c_e_combined)
+        np.testing.assert_array_almost_equal(self.c_e(t, x), c_e_combined, decimal=14)
 
     def test_all(self):
         self.test_concentration_limit()

From c2c0f28e6d23691adf73b33dc46b784235b694b8 Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Wed, 9 Nov 2022 16:17:39 -0500
Subject: [PATCH 11/12] coverage

---
 pybamm/expression_tree/concatenations.py                | 2 --
 pybamm/meshes/meshes.py                                 | 3 ---
 tests/unit/test_expression_tree/test_unary_operators.py | 8 ++++++++
 3 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 01b2ea870b..e575dc3110 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -437,8 +437,6 @@ def simplified_concatenation(*children):
 def concatenation(*children):
     """Helper function to create concatenations."""
     # TODO: add option to turn off simplifications
-    if len(children) == 3 and children[-1].domain == ["current collector"]:
-        print("here")
     return simplified_concatenation(*children)
 
 
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 5e1bb0cd7a..793ceced44 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -152,9 +152,6 @@ def combine_submeshes(self, *submeshnames):
         """
         if submeshnames == ():
             raise ValueError("Submesh domains being combined cannot be empty")
-        # If there is just a single submesh, we can return it directly
-        if len(submeshnames) == 1:
-            return self[submeshnames[0]]
         # Check that the final edge of each submesh is the same as the first edge of the
         # next submesh
         for i in range(len(submeshnames) - 1):
diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py
index 9dd36d4d7d..7d4906558f 100644
--- a/tests/unit/test_expression_tree/test_unary_operators.py
+++ b/tests/unit/test_expression_tree/test_unary_operators.py
@@ -651,6 +651,14 @@ def test_to_equation(self):
             sympy.Integral("d", sympy.Symbol("xn")),
         )
 
+    def test_explicit_time_integral(self):
+        expr = pybamm.ExplicitTimeIntegral(pybamm.Parameter("param"), pybamm.Scalar(1))
+        self.assertEqual(expr.child, pybamm.Parameter("param"))
+        self.assertEqual(expr.initial_condition, pybamm.Scalar(1))
+        self.assertEqual(expr.name, "explicit time integral")
+        self.assertEqual(expr.new_copy(), expr)
+        self.assertFalse(expr.is_constant())
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

From 37b96d02d1e1753d40a46b052a3291b017f96a7e Mon Sep 17 00:00:00 2001
From: Valentin Sulzer <valentinsulzer@hotmail.com>
Date: Thu, 10 Nov 2022 08:16:13 -0500
Subject: [PATCH 12/12] Update CHANGELOG.md

---
 CHANGELOG.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3926682e27..f32df3e538 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,7 +6,7 @@
 
 ## Optimizations
 
-- Added more rules for simplifying expressions, especially around Concatenations. Also, meshes constructed from multiple domains are now memoized ([#2443](https://github.com/pybamm-team/PyBaMM/pull/2443))
+- Added more rules for simplifying expressions, especially around Concatenations. Also, meshes constructed from multiple domains are now cached ([#2443](https://github.com/pybamm-team/PyBaMM/pull/2443))
 - Added more rules for simplifying expressions. Constants in binary operators are now moved to the left by default (e.g. `x*2` returns `2*x`) ([#2424](https://github.com/pybamm-team/PyBaMM/pull/2424))
 
 ## Breaking changes