Repository: tinosulzer/PyBaMM

Issue:
One Particle vs Many Particles
We need a method to distinguish the cases of having 1 particle in an electrode and many particles in an electrode. This will also be relevant for when we have 2D current collectors and 1D DFN models between.

**Possible Approach**

1) State the spatial variables and their limits (this step just makes it easier to put the same spatial variables into the geometry): 
```spatial_vars = {xn: {"min": 0, "max": ln}, .... etc} ``` 

2) Assign spatial variables as either primary or secondary to each domain
```domain_vars = {"negative electrode": {"primary": [xn, y, z], "secondary": []}, "negative particle": {"primary": [r], "secondary": [xn]}}```

3) Process into same form as before except with primary and secondary variables
```geometry = Geometry(domain_vars, spatial_vars) ```
form: 
```geometry = {"negative electrode": {"primary": {x: {"min": 0, "max": ln}, "secondary": []}, "negative particle": {"primary": {rn: {"min": 0, "max": 1} , "secondary": {xn: {"min": 0, "max": ln"} )]}}```

4) 2 minor changes to `Mesh`
    - scan through both primary and secondary variables and create `SubMesh`
    - create a flag which is a list of length equal to the dimension of the `SubMesh` matrix with entries `["primary", "primary" "secondary"]` corresponding to the type of variable each dimension is associated with. 

5) Change to `Discretisation`: 
    - scan the `SubMesh` for the `["primary", "primary" "secondary"]` flags and also assign them to the spatial methods upon initialisation
 
6) Change to `SpatialMethod`: 
    - change grad/div so that they read the assigned  `["primary", "primary" "secondary"]` flag and apply the derivatives in the correct way. 

This method is general so can be applied to both the many particle problem and the case where we have 1D DFN models coupled to 2D current collectors. 

**Downside** 
To some extent we are hiding what may be regarded as a part of the model equations (that grad should just be d/dr) within the Geometry. But, if we don't do this we would have to add more symbols for d/dr etc and also re-write the same equations for different dimensions. 

 **Motivation**
Required to implement the SPMe, DFN, and DFNCC(2+1+1D). 


Discussion:

Patch:
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 6f87435a..464d9cb5 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -132,6 +132,7 @@ from .geometry.geometry import (
     Geometry,
     Geometry1DMacro,
     Geometry1DMicro,
+    Geometry1p1DMicro,
     Geometry3DMacro,
 )
 
diff --git a/pybamm/discretisations/discretisation.py b/pybamm/discretisations/discretisation.py
index 0e132481..b20b3694 100644
--- a/pybamm/discretisations/discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -103,7 +103,8 @@ class Discretisation(object):
             # Otherwise, add up the size of all the domains in variable.domain
             else:
                 for dom in variable.domain:
-                    end += self._spatial_methods[dom].mesh[dom].npts_for_broadcast
+                    for submesh in self._spatial_methods[dom].mesh[dom]:
+                        end += submesh.npts_for_broadcast
             y_slices[variable.id] = slice(start, end)
             start = end
         self._y_slices = y_slices
@@ -218,7 +219,7 @@ class Discretisation(object):
             child = symbol.children[0]
             discretised_child = self.process_symbol(child)
             return self._spatial_methods[child.domain[0]].integral(
-                child.domain, discretised_child
+                child.domain, child, discretised_child
             )
 
         elif isinstance(symbol, pybamm.Broadcast):
diff --git a/pybamm/expression_tree/binary_operators.py b/pybamm/expression_tree/binary_operators.py
index 432d0088..27977202 100644
--- a/pybamm/expression_tree/binary_operators.py
+++ b/pybamm/expression_tree/binary_operators.py
@@ -233,6 +233,7 @@ class MatrixMultiplication(BinaryOperator):
 
     def evaluate(self, t=None, y=None):
         """ See :meth:`pybamm.Symbol.evaluate()`. """
+
         return self.children[0].evaluate(t, y) @ self.children[1].evaluate(t, y)
 
     def _binary_simplify(self, left, right):
diff --git a/pybamm/expression_tree/broadcasts.py b/pybamm/expression_tree/broadcasts.py
index ec98e261..5397f6df 100644
--- a/pybamm/expression_tree/broadcasts.py
+++ b/pybamm/expression_tree/broadcasts.py
@@ -76,9 +76,12 @@ class NumpyBroadcast(Broadcast):
         if domain == []:
             self.broadcasting_vector_size = 1
         else:
-            self.broadcasting_vector_size = sum(
-                [mesh[dom].npts_for_broadcast for dom in domain]
-            )
+            vector_size = 0
+            for dom in domain:
+                # just create a vector of the points even in 2 and 3D
+                for i in range(len(mesh[dom])):
+                    vector_size += mesh[dom][i].npts_for_broadcast
+            self.broadcasting_vector_size = vector_size
         # create broadcasting vector (vector of ones with shape determined by the
         # domain)
         self.broadcasting_vector = np.ones(self.broadcasting_vector_size)
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index ab2ddf7e..a6dcf5ab 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -137,7 +137,9 @@ class DomainConcatenation(Concatenation):
         start = 0
         end = 0
         for dom in node.domain:
-            end += mesh[dom].npts
+            prim_pts = mesh[dom][0].npts
+            second_pts = len(mesh[dom])
+            end += prim_pts * second_pts
             slices[dom] = slice(start, end)
             start = end
         return slices
diff --git a/pybamm/geometry/geometry.py b/pybamm/geometry/geometry.py
index ccf007bf..e4828a4c 100644
--- a/pybamm/geometry/geometry.py
+++ b/pybamm/geometry/geometry.py
@@ -8,7 +8,14 @@ import pybamm
 
 
 class Geometry(dict):
-    """A geometry class to store the details features of the cell geometry
+    """A geometry class to store the details features of the cell geometry.
+        Geometry extends the class dictionary and uses the key words:
+        "negative electrode", "positive electrode", etc to indicate the subdomain.
+        Within each subdomain, there are "primary" and "secondary" dimensions.
+        "primary" dimensions correspond to dimensiones on which spatial
+        operators will be applied (e.g. the gradient and divergence). In contrast,
+        spatial operators do not act along "secondary" dimensions. This allows for
+        multiple independent particles to be included into a model.
 
      **Extends**: :class:`dict`
 
@@ -26,6 +33,8 @@ class Geometry(dict):
                 geometry = Geometry3DMacro()
             elif geometry == "1D micro":
                 geometry = Geometry1DMicro()
+            elif geometry == "1+1D micro":
+                geometry = Geometry1p1DMicro()
             # avoid combining geometries that clash
             if any([k in self.keys() for k in geometry.keys()]):
                 raise ValueError("trying to overwrite existing geometry")
@@ -42,9 +51,13 @@ class Geometry1DMacro(Geometry):
         l_n = pybamm.standard_parameters.l_n
         l_s = pybamm.standard_parameters.l_s
 
-        self["negative electrode"] = {x: {"min": pybamm.Scalar(0), "max": l_n}}
-        self["separator"] = {x: {"min": l_n, "max": l_n + l_s}}
-        self["positive electrode"] = {x: {"min": l_n + l_s, "max": pybamm.Scalar(1)}}
+        self["negative electrode"] = {
+            "primary": {x: {"min": pybamm.Scalar(0), "max": l_n}}
+        }
+        self["separator"] = {"primary": {x: {"min": l_n, "max": l_n + l_s}}}
+        self["positive electrode"] = {
+            "primary": {x: {"min": l_n + l_s, "max": pybamm.Scalar(1)}}
+        }
 
         # update with custom geometry if non empty
         self.update(custom_geometry)
@@ -58,10 +71,33 @@ class Geometry1DMicro(Geometry):
         r_p = pybamm.SpatialVariable("r", ["positive particle"])
 
         self["negative particle"] = {
-            r_n: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
+            "primary": {r_n: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}
+        }
+        self["positive particle"] = {
+            "primary": {r_p: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}
+        }
+        # update with custom geometry if non empty
+        self.update(custom_geometry)
+
+
+class Geometry1p1DMicro(Geometry):
+    def __init__(self, custom_geometry={}):
+        super().__init__()
+
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        x = pybamm.SpatialVariable("x", whole_cell)
+        r_n = pybamm.SpatialVariable("r", ["negative particle"])
+        r_p = pybamm.SpatialVariable("r", ["positive particle"])
+        l_n = pybamm.standard_parameters.l_n
+        l_s = pybamm.standard_parameters.l_s
+
+        self["negative particle"] = {
+            "primary": {r_n: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}},
+            "secondary": {x: {"min": pybamm.Scalar(0), "max": l_n}},
         }
         self["positive particle"] = {
-            r_p: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
+            "primary": {r_p: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}},
+            "secondary": {x: {"min": l_n + l_s, "max": pybamm.Scalar(1)}},
         }
         # update with custom geometry if non empty
         self.update(custom_geometry)
@@ -91,6 +127,6 @@ class Geometry3DMacro(Geometry1DMacro):
 
         MACRO_DOMAINS = ["negative electrode", "separator", "positive electrode"]
         for domain in MACRO_DOMAINS:
-            self[domain][y] = y_lim
-            self[domain][z] = z_lim
+            self[domain]["primary"][y] = y_lim
+            self[domain]["primary"][z] = z_lim
         self.update(custom_geometry)
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index 14392cdf..64909a1d 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -19,7 +19,7 @@ KNOWN_DOMAINS = [
 
 class Mesh(dict):
     """
-    Mesh contains the submeshes
+    Mesh contains a list of submeshes on each subdomain.
 
     **Extends**: dict
 
@@ -27,7 +27,7 @@ class Mesh(dict):
     ----------
 
     geometry : :class: `Geometry`
-        contains the geometry of the problem
+        contains the geometry of the problem.
     submesh_types: dict
         contains the types of submeshes to use (e.g. Uniform1DSubMesh)
     submesh_pts: dict
@@ -39,9 +39,13 @@ class Mesh(dict):
         super().__init__()
         self.submesh_pts = submesh_pts
         for domain in geometry:
-            submesh_type = submesh_types[domain]
-            submesh_pt = submesh_pts[domain]
-            self[domain] = submesh_type(geometry[domain], submesh_pt)
+            repeats = 1
+            if "secondary" in geometry[domain].keys():
+                for var in geometry[domain]["secondary"].keys():
+                    repeats = submesh_pts[domain][var.name]  # note (specific to FV)
+            self[domain] = [
+                submesh_types[domain](geometry[domain]["primary"], submesh_pts[domain])
+            ] * repeats
 
     def combine_submeshes(self, *submeshnames):
         """Combine submeshes into a new submesh, using self.submeshclass
@@ -60,24 +64,22 @@ class Mesh(dict):
         """
         # Check that the final edge of each submesh is the same as the first edge of the
         # next submesh
-
-        edges_aligned = all(
-            [
-                self[submeshnames[i]].edges[-1] == self[submeshnames[i + 1]].edges[0]
-                for i in range(len(submeshnames) - 1)
-            ]
-        )
-
-        if edges_aligned:
-            # Combine submeshes, being careful not to double-count repeated edges at the
-            # intersection of submeshes
+        for i in range(len(submeshnames) - 1):
+            for j in range(len(self[submeshnames[i]])):
+                if (
+                    self[submeshnames[i]][j].edges[-1]
+                    != self[submeshnames[i + 1]][j].edges[0]
+                ):
+                    raise pybamm.DomainError("submesh edges are not aligned")
+
+        submeshes = [None] * len(self[submeshnames[0]])
+        for i in range(len(self[submeshnames[0]])):
             combined_submesh_edges = np.concatenate(
-                [self[submeshnames[0]].edges]
-                + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]]
+                [self[submeshnames[0]][i].edges]
+                + [self[submeshname][i].edges[1:] for submeshname in submeshnames[1:]]
             )
-            return pybamm.SubMesh1D(combined_submesh_edges)
-        else:
-            raise pybamm.DomainError("submesh edges are not aligned")
+            submeshes[i] = pybamm.SubMesh1D(combined_submesh_edges)
+        return submeshes
 
     def add_ghost_meshes(self):
         """
@@ -87,15 +89,22 @@ class Mesh(dict):
         """
         # Get all submeshes relating to space (i.e. exclude time)
         submeshes = [
-            (name, submesh) for name, submesh in self.items() if name != "time"
+            (domain, submesh_list)
+            for domain, submesh_list in self.items()
+            if domain != "time"
         ]
-        for submeshname, submesh in submeshes:
-            edges = submesh.edges
+        for domain, submesh_list in submeshes:
+
+            self[domain + "_left ghost cell"] = [None] * len(submesh_list)
+            self[domain + "_right ghost cell"] = [None] * len(submesh_list)
+            for i, submesh in enumerate(submesh_list):
+                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[submeshname + "_left ghost cell"] = pybamm.SubMesh1D(lgs_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"][i] = pybamm.SubMesh1D(lgs_edges)
 
-            # 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[submeshname + "_right ghost cell"] = pybamm.SubMesh1D(rgs_edges)
+                # 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"][i] = pybamm.SubMesh1D(rgs_edges)
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index 8f27f5be..e00c226e 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -125,11 +125,12 @@ class ParameterValues(dict):
             """
 
         for domain in geometry:
-            for spatial_variable, spatial_limits in geometry[domain].items():
-                for lim, sym in spatial_limits.items():
-                    geometry[domain][spatial_variable][lim] = self.process_symbol(
-                        sym
-                    ).evaluate()
+            for prim_sec, variables in geometry[domain].items():
+                for spatial_variable, spatial_limits in variables.items():
+                    for lim, sym in spatial_limits.items():
+                        geometry[domain][prim_sec][spatial_variable][
+                            lim
+                        ] = self.process_symbol(sym).evaluate()
 
     def process_symbol(self, symbol):
         """Walk through the symbol and replace any Parameter with a Value.
diff --git a/pybamm/spatial_methods/finite_volume.py b/pybamm/spatial_methods/finite_volume.py
index b6ebea15..6714f70f 100644
--- a/pybamm/spatial_methods/finite_volume.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -7,6 +7,8 @@ import pybamm
 
 import numpy as np
 from scipy.sparse import spdiags
+from scipy.sparse import eye
+from scipy.sparse import kron
 
 
 class FiniteVolume(pybamm.SpatialMethod):
@@ -25,7 +27,8 @@ class FiniteVolume(pybamm.SpatialMethod):
     def __init__(self, mesh):
         # add npts_for_broadcast to mesh domains for this particular discretisation
         for dom in mesh.keys():
-            mesh[dom].npts_for_broadcast = mesh[dom].npts
+            for i in range(len(mesh[dom])):
+                mesh[dom][i].npts_for_broadcast = mesh[dom][i].npts
         super().__init__(mesh)
 
     def spatial_variable(self, symbol):
@@ -46,7 +49,7 @@ class FiniteVolume(pybamm.SpatialMethod):
         # for finite volume we use the cell centres
         if symbol.name in ["x", "r"]:
             symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
-            return pybamm.Vector(symbol_mesh.nodes)
+            return pybamm.Vector(symbol_mesh[0].nodes)
         else:
             raise NotImplementedError("3D meshes not yet implemented")
 
@@ -58,7 +61,6 @@ class FiniteVolume(pybamm.SpatialMethod):
         See :meth: `pybamm.SpatialMethod.broadcast`
         """
 
-        # for finite volume we send variables to cells and so use number_of_cells
         broadcasted_symbol = pybamm.NumpyBroadcast(symbol, domain, self.mesh)
 
         # if the broadcasted symbol evaluates to a constant value, replace the
@@ -85,7 +87,9 @@ class FiniteVolume(pybamm.SpatialMethod):
         if symbol.id in boundary_conditions:
             lbc = boundary_conditions[symbol.id]["left"]
             rbc = boundary_conditions[symbol.id]["right"]
-            discretised_symbol = self.add_ghost_nodes(discretised_symbol, lbc, rbc)
+            discretised_symbol = self.add_ghost_nodes(
+                symbol, discretised_symbol, lbc, rbc
+            )
             domain = (
                 [domain[0] + "_left ghost cell"]
                 + domain
@@ -112,16 +116,26 @@ class FiniteVolume(pybamm.SpatialMethod):
             The (sparse) finite volume gradient matrix for the domain
         """
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh_list = self.mesh.combine_submeshes(*domain)
 
-        # Create matrix using submesh
+        # can just use 1st entry of list to obtain the point etc
+        submesh = submesh_list[0]
+
+        # Create 1D matrix using submesh
         n = submesh.npts
         e = 1 / submesh.d_nodes
         data = np.vstack(
             [np.concatenate([-e, np.array([0])]), np.concatenate([np.array([0]), e])]
         )
         diags = np.array([0, 1])
-        matrix = spdiags(data, diags, n - 1, n)
+        sub_matrix = spdiags(data, diags, n - 1, n)
+
+        # second dim length
+        second_dim_len = len(submesh_list)
+
+        # generate full matrix from the submatrix
+        matrix = kron(eye(second_dim_len), sub_matrix)
+
         return pybamm.Matrix(matrix)
 
     def divergence(self, symbol, discretised_symbol, boundary_conditions):
@@ -133,30 +147,70 @@ class FiniteVolume(pybamm.SpatialMethod):
             assert isinstance(key, int), TypeError(
                 "boundary condition keys should be hashes, not {}".format(type(key))
             )
+
+        domain = symbol.domain
+        submesh_list = self.mesh.combine_submeshes(*domain)
+        # create a bc vector of length equal to the number variables
+        # (only has non zero entries for neumann bcs)
+
+        prim_dim = submesh_list[0].npts
+        second_dim = len(submesh_list)
+        total_pts = prim_dim * second_dim
+
         # Add Neumann boundary conditions if defined
         if symbol.id in boundary_conditions:
-            # for the particles there will be a "negative particle" "left" and "right"
-            # and also a "positive particle" left and right.
+            # TODO:these are symbols so need to check them
             lbc = boundary_conditions[symbol.id]["left"]
             rbc = boundary_conditions[symbol.id]["right"]
-            discretised_symbol = pybamm.NumpyConcatenation(lbc, discretised_symbol, rbc)
 
-        domain = symbol.domain
+            # doing via loop so that it is easier to implement x varing bcs
+            bcs_symbol = pybamm.Vector(np.array([]))  # empty vector
+            for i in range(len(submesh_list)):
+                # only the interior equations:
+                interior = pybamm.Vector(np.zeros(prim_dim - 2))
+                left = -lbc / pybamm.Vector(np.array([submesh_list[i].d_edges[0]]))
+                right = rbc / pybamm.Vector(np.array([submesh_list[i].d_edges[-1]]))
+                bcs_symbol = pybamm.NumpyConcatenation(
+                    bcs_symbol, left, interior, right
+                )
+
+            # now we must create a matrix of size (npts * (npts -1) )
+            # this is a different size to the one created when we have
+            # flux boundary conditions so need a flag
+            divergence_matrix = self.divergence_matrix(domain, bc_type="neumann")
+
+            # only need interior edges for spherical neumann
+            edges = submesh_list[0].edges[1:-1]
+
+        else:
+            divergence_matrix = self.divergence_matrix(domain)
+            bcs_vec = np.zeros(total_pts)
+            bcs_symbol = pybamm.Vector(bcs_vec)
+            # need all edges for spherical dirichlet
+            edges = submesh_list[0].edges
+
         # check for particle domain
-        divergence_matrix = self.divergence_matrix(domain)
         if ("negative particle" or "positive particle") in domain:
-            submesh = self.mesh.combine_submeshes(*domain)
-            r = pybamm.Vector(submesh.nodes)
-            r_edges = pybamm.Vector(submesh.edges)
+
+            # create np.array of repeated submesh[0].nodes
+            r_numpy = np.kron(np.ones(second_dim), submesh_list[0].nodes)
+            r_edges_numpy = np.kron(np.ones(second_dim), edges)
+
+            r = pybamm.Vector(r_numpy)
+            r_edges = pybamm.Vector(r_edges_numpy)
+
+            # for clarity, we are implicitly multiplying the the lbc by r^2=0
+            # and the rbc by r^2=1. But lbc is 0 so we don't need to do
+            # any r_edges^2 operations on bcs_symbol
             out = (1 / (r ** 2)) * (
-                divergence_matrix @ ((r_edges ** 2) * discretised_symbol)
+                divergence_matrix @ ((r_edges ** 2) * discretised_symbol) + bcs_symbol
             )
         else:
-            divergence_matrix = self.divergence_matrix(domain)
-            out = divergence_matrix @ discretised_symbol
+            out = divergence_matrix @ discretised_symbol + bcs_symbol
+
         return out
 
-    def divergence_matrix(self, domain):
+    def divergence_matrix(self, domain, bc_type="dirichlet"):
         """
         Divergence matrix for finite volumes in the appropriate domain.
         Equivalent to div(N) = (N[1:] - N[:-1])/dx
@@ -172,28 +226,52 @@ class FiniteVolume(pybamm.SpatialMethod):
             The (sparse) finite volume divergence matrix for the domain
         """
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh_list = self.mesh.combine_submeshes(*domain)
+
+        # can just use 1st entry of list to obtain the point etc
+        submesh = submesh_list[0]
+        e = 1 / submesh.d_edges
 
         # Create matrix using submesh
         n = submesh.npts + 1
-        e = 1 / submesh.d_edges
-        data = np.vstack(
-            [np.concatenate([-e, np.array([0])]), np.concatenate([np.array([0]), e])]
-        )
-        diags = np.array([0, 1])
-        matrix = spdiags(data, diags, n - 1, n)
+        if bc_type == "dirichlet":
+            data = np.vstack(
+                [
+                    np.concatenate([-e, np.array([0])]),
+                    np.concatenate([np.array([0]), e]),
+                ]
+            )
+            diags = np.array([0, 1])
+            sub_matrix = spdiags(data, diags, n - 1, n)
+        elif bc_type == "neumann":
+            # we don't have to act on bc fluxes which are now in
+            # the bc vector
+            data = np.vstack([-e[1:], e[:-1]])
+            diags = np.array([-1, 0])
+            sub_matrix = spdiags(data, diags, n - 1, n - 2)
+        else:
+            raise NotImplementedError(
+                "Can only process Neumann or Dirichlet boundary conditions"
+            )
+
+        # repeat matrix for each node in secondary dimensions
+        second_dim_len = len(submesh_list)
+        # generate full matrix from the submatrix
+        matrix = kron(eye(second_dim_len), sub_matrix)
         return pybamm.Matrix(matrix)
 
-    def integral(self, domain, discretised_symbol):
+    def integral(self, domain, symbol, discretised_symbol):
         """Vector-vector dot product to implement the integral operator.
         See :meth:`pybamm.BaseDiscretisation.integral`
         """
         # Calculate integration vector
         integration_vector = self.definite_integral_vector(domain)
         # Check for particle domain
-        if ("negative particle" or "positive particle") in domain:
-            submesh = self.mesh.combine_submeshes(*domain)
-            r = pybamm.Vector(submesh.nodes)
+        if ("negative particle" or "positive particle") in symbol.domain:
+            submesh_list = self.mesh.combine_submeshes(*symbol.domain)
+            second_dim = len(submesh_list)
+            r_numpy = np.kron(np.ones(second_dim), submesh_list[0].nodes)
+            r = pybamm.Vector(r_numpy)
             out = 2 * np.pi * integration_vector @ (discretised_symbol * r)
         else:
             out = integration_vector @ discretised_symbol
@@ -221,14 +299,16 @@ class FiniteVolume(pybamm.SpatialMethod):
             The finite volume integral vector for the domain
         """
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh_list = self.mesh.combine_submeshes(*domain)
 
         # Create vector of ones using submesh
-        vector = submesh.d_edges * np.ones_like(submesh.nodes)
+        vector = np.array([])
+        for submesh in submesh_list:
+            vector = np.append(vector, submesh.d_edges * np.ones_like(submesh.nodes))
 
         return pybamm.Vector(vector)
 
-    def add_ghost_nodes(self, discretised_symbol, lbc, rbc):
+    def add_ghost_nodes(self, symbol, discretised_symbol, lbc, rbc):
         """
         Add Dirichlet boundary conditions via ghost nodes.
 
@@ -265,18 +345,47 @@ class FiniteVolume(pybamm.SpatialMethod):
                 type(discretised_symbol)
             )
         )
-        # left ghost cell
+
+        # determine the y_slice sizes
         y_slice_start = discretised_symbol.y_slice.start
-        first_node = pybamm.StateVector(slice(y_slice_start, y_slice_start + 1))
-        left_ghost_cell = 2 * lbc - first_node
-        # right ghost cell
         y_slice_stop = discretised_symbol.y_slice.stop
-        last_node = pybamm.StateVector(slice(y_slice_stop - 1, y_slice_stop))
-        right_ghost_cell = 2 * rbc - last_node
-        # concatenate
-        return pybamm.NumpyConcatenation(
-            left_ghost_cell, discretised_symbol, right_ghost_cell
-        )
+
+        y = np.arange(y_slice_start, y_slice_stop)
+
+        # reshape y_slices into more helpful form
+        submesh_list = self.mesh.combine_submeshes(*symbol.domain)
+        if isinstance(submesh_list[0].npts, list):
+            NotImplementedError("Can only take in 1D primary directions")
+        # size = [submesh_list[0].npts, len(submesh_list)]
+        size = [len(submesh_list), submesh_list[0].npts]
+        y = np.reshape(y, size)
+        y_left = y[:, 0]
+        y_right = y[:, -1]
+
+        new_discretised_symbol = pybamm.Vector(np.array([]))  # starts empty
+
+        for i in range(len(submesh_list)):
+            y_slice_start = y_left[i]
+            y_slice_stop = y_right[i]
+
+            # left ghost cell
+            first_node = pybamm.StateVector(slice(y_slice_start, y_slice_start + 1))
+            left_ghost_cell = 2 * lbc - first_node
+            # middle symbol
+            sub_disc_symbol = pybamm.StateVector(slice(y_slice_start, y_slice_stop + 1))
+            # right ghost cell
+            last_node = pybamm.StateVector(slice(y_slice_stop, y_slice_stop + 1))
+            right_ghost_cell = 2 * rbc - last_node
+
+            concatenated_sub_disc_symbol = pybamm.NumpyConcatenation(
+                left_ghost_cell, sub_disc_symbol, right_ghost_cell
+            )
+
+            new_discretised_symbol = pybamm.NumpyConcatenation(
+                new_discretised_symbol, concatenated_sub_disc_symbol
+            )
+
+        return new_discretised_symbol
 
     def surface_value(self, discretised_symbol):
         """


Test patch:
diff --git a/tests/__init__.py b/tests/__init__.py
index 37653849..1d15c629 100644
--- a/tests/__init__.py
+++ b/tests/__init__.py
@@ -6,4 +6,8 @@ from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
 from .test_models.standard_model_tests import StandardModelTest
-from .shared import get_mesh_for_testing, get_discretisation_for_testing
+from .shared import (
+    get_mesh_for_testing,
+    get_p2d_mesh_for_testing,
+    get_discretisation_for_testing,
+)
diff --git a/tests/shared.py b/tests/shared.py
index b13ab334..ac4206e7 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -11,13 +11,14 @@ class SpatialMethodForTesting(pybamm.SpatialMethod):
 
     def __init__(self, mesh):
         for dom in mesh.keys():
-            mesh[dom].npts_for_broadcast = mesh[dom].npts
+            for i in range(len(mesh[dom])):
+                mesh[dom][i].npts_for_broadcast = mesh[dom][i].npts
         super().__init__(mesh)
 
     def spatial_variable(self, symbol):
         # for finite volume we use the cell centres
         symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
-        return pybamm.Vector(symbol_mesh.nodes)
+        return pybamm.Vector(symbol_mesh[0].nodes)
 
     def broadcast(self, symbol, domain):
         # for finite volume we send variables to cells and so use number_of_cells
@@ -35,14 +36,14 @@ class SpatialMethodForTesting(pybamm.SpatialMethod):
     def gradient(self, symbol, discretised_symbol, boundary_conditions):
         n = 0
         for domain in symbol.domain:
-            n += self.mesh[domain].npts
+            n += self.mesh[domain][0].npts
         gradient_matrix = pybamm.Matrix(np.eye(n))
         return gradient_matrix @ discretised_symbol
 
     def divergence(self, symbol, discretised_symbol, boundary_conditions):
         n = 0
         for domain in symbol.domain:
-            n += self.mesh[domain].npts
+            n += self.mesh[domain][0].npts
         divergence_matrix = pybamm.Matrix(np.eye(n))
         return divergence_matrix @ discretised_symbol
 
@@ -90,6 +91,48 @@ def get_mesh_for_testing(npts=None):
     return pybamm.Mesh(geometry, submesh_types, submesh_pts)
 
 
+def get_p2d_mesh_for_testing(npts=None, mpts=None):
+    param = pybamm.ParameterValues(
+        base_parameters={
+            "Negative electrode width": 0.3,
+            "Separator width": 0.2,
+            "Positive electrode width": 0.3,
+        }
+    )
+
+    geometry = pybamm.Geometry("1D macro", "1+1D micro")
+    param.process_geometry(geometry)
+
+    # provide mesh properties
+    submesh_types = {
+        "negative electrode": pybamm.Uniform1DSubMesh,
+        "separator": pybamm.Uniform1DSubMesh,
+        "positive electrode": pybamm.Uniform1DSubMesh,
+        "negative particle": pybamm.Uniform1DSubMesh,
+        "positive particle": pybamm.Uniform1DSubMesh,
+    }
+
+    if mpts is None:
+        submesh_pts = {
+            "negative electrode": {"x": 40},
+            "separator": {"x": 25},
+            "positive electrode": {"x": 35},
+            "negative particle": {"r": 10, "x": 40},
+            "positive particle": {"r": 10, "x": 35},
+        }
+    else:
+        n = 3 * round(npts / 3)
+        submesh_pts = {
+            "negative electrode": {"x": n},
+            "separator": {"x": n},
+            "positive electrode": {"x": n},
+            "negative particle": {"r": mpts, "x": n},
+            "positive particle": {"r": mpts, "x": n},
+        }
+
+    return pybamm.Mesh(geometry, submesh_types, submesh_pts)
+
+
 def get_discretisation_for_testing(npts=None):
     mesh = get_mesh_for_testing(npts)
     spatial_methods = {
diff --git a/tests/test_discretisations/test_discretisation.py b/tests/test_discretisations/test_discretisation.py
index bdd75122..418a4af1 100644
--- a/tests/test_discretisations/test_discretisation.py
+++ b/tests/test_discretisations/test_discretisation.py
@@ -50,7 +50,7 @@ class TestDiscretise(unittest.TestCase):
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
-        c_true = combined_submesh.nodes ** 2
+        c_true = combined_submesh[0].nodes ** 2
         y = c_true
         np.testing.assert_array_equal(y[disc._y_slices[c.id]], c_true)
 
@@ -64,8 +64,8 @@ class TestDiscretise(unittest.TestCase):
             disc._y_slices,
             {c.id: slice(0, 100), d.id: slice(100, 200), jn.id: slice(200, 240)},
         )
-        d_true = 4 * combined_submesh.nodes
-        jn_true = mesh["negative electrode"].nodes ** 3
+        d_true = 4 * combined_submesh[0].nodes
+        jn_true = mesh["negative electrode"][0].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
         np.testing.assert_array_equal(y[disc._y_slices[c.id]], c_true)
         np.testing.assert_array_equal(y[disc._y_slices[d.id]], d_true)
@@ -87,8 +87,8 @@ class TestDiscretise(unittest.TestCase):
                 jp.id: slice(265, 300),
             },
         )
-        d_true = 4 * combined_submesh.nodes
-        jn_true = mesh["negative electrode"].nodes ** 3
+        d_true = 4 * combined_submesh[0].nodes
+        jn_true = mesh["negative electrode"][0].nodes ** 3
         y = np.concatenate([c_true, d_true, jn_true])
         np.testing.assert_array_equal(y[disc._y_slices[c.id]], c_true)
         np.testing.assert_array_equal(y[disc._y_slices[d.id]], d_true)
@@ -210,7 +210,7 @@ class TestDiscretise(unittest.TestCase):
             self.assertIsInstance(eqn_disc.children[1], pybamm.StateVector)
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
-            y = combined_submesh.nodes ** 2
+            y = combined_submesh[0].nodes ** 2
             var_disc = disc.process_symbol(var)
             # grad and var are identity operators here (for testing purposes)
             np.testing.assert_array_equal(
@@ -227,7 +227,7 @@ class TestDiscretise(unittest.TestCase):
             self.assertIsInstance(eqn_disc.children[1].children[0], pybamm.Matrix)
             self.assertIsInstance(eqn_disc.children[1].children[1], pybamm.StateVector)
 
-            y = combined_submesh.nodes ** 2
+            y = combined_submesh[0].nodes ** 2
             var_disc = disc.process_symbol(var)
             # grad and var are identity operators here (for testing purposes)
             np.testing.assert_array_equal(
@@ -263,7 +263,7 @@ class TestDiscretise(unittest.TestCase):
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
-        y = combined_submesh.nodes ** 2
+        y = combined_submesh[0].nodes ** 2
         disc._bcs = boundary_conditions
 
         disc.set_variable_slices(list(rhs.keys()))
@@ -273,7 +273,7 @@ class TestDiscretise(unittest.TestCase):
         # initial conditions
         y0 = disc.process_dict(initial_conditions)
         np.testing.assert_array_equal(
-            y0[c].evaluate(0, None), 3 * np.ones_like(combined_submesh.nodes)
+            y0[c].evaluate(0, None), 3 * np.ones_like(combined_submesh[0].nodes)
         )
         # vars
         processed_vars = disc.process_dict(variables)
@@ -288,7 +288,7 @@ class TestDiscretise(unittest.TestCase):
         initial_conditions = {c: pybamm.Scalar(3), T: pybamm.Scalar(5)}
         boundary_conditions = {}
         y = np.concatenate(
-            [combined_submesh.nodes ** 2, mesh["negative electrode"].nodes ** 4]
+            [combined_submesh[0].nodes ** 2, mesh["negative electrode"][0].nodes ** 4]
         )
 
         variables = list(rhs.keys())
@@ -304,10 +304,11 @@ class TestDiscretise(unittest.TestCase):
         # initial conditions
         y0 = disc.process_dict(initial_conditions)
         np.testing.assert_array_equal(
-            y0[c].evaluate(0, None), 3 * np.ones_like(combined_submesh.nodes)
+            y0[c].evaluate(0, None), 3 * np.ones_like(combined_submesh[0].nodes)
         )
         np.testing.assert_array_equal(
-            y0[T].evaluate(0, None), 5 * np.ones_like(mesh["negative electrode"].nodes)
+            y0[T].evaluate(0, None),
+            5 * np.ones_like(mesh["negative electrode"][0].nodes),
         )
 
     def test_process_variables_dict(self):
@@ -340,7 +341,7 @@ class TestDiscretise(unittest.TestCase):
         disc.process_model(model)
 
         y0 = model.concatenated_initial_conditions
-        np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh.nodes))
+        np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh[0].nodes))
         np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
         # grad and div are identity operators here
         np.testing.assert_array_equal(y0, model.variables["c"].evaluate(None, y0))
@@ -371,16 +372,17 @@ class TestDiscretise(unittest.TestCase):
             y0,
             np.concatenate(
                 [
-                    2 * np.ones_like(combined_submesh.nodes),
-                    5 * np.ones_like(mesh["negative electrode"].nodes),
-                    8 * np.ones_like(mesh["negative electrode"].nodes),
+                    2 * np.ones_like(combined_submesh[0].nodes),
+                    5 * np.ones_like(mesh["negative electrode"][0].nodes),
+                    8 * np.ones_like(mesh["negative electrode"][0].nodes),
                 ]
             ),
         )
         # grad and div are identity operators here
         np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
         c0, T0, S0 = np.split(
-            y0, np.cumsum([combined_submesh.npts, mesh["negative electrode"].npts])
+            y0,
+            np.cumsum([combined_submesh[0].npts, mesh["negative electrode"][0].npts]),
         )
         np.testing.assert_array_equal(S0 * T0, model.variables["ST"].evaluate(None, y0))
 
@@ -421,20 +423,20 @@ class TestDiscretise(unittest.TestCase):
             y0,
             np.concatenate(
                 [
-                    3 * np.ones_like(combined_submesh.nodes),
-                    6 * np.ones_like(combined_submesh.nodes),
+                    3 * np.ones_like(combined_submesh[0].nodes),
+                    6 * np.ones_like(combined_submesh[0].nodes),
                 ]
             ),
         )
 
         # grad and div are identity operators here
         np.testing.assert_array_equal(
-            y0[: combined_submesh.npts], model.concatenated_rhs.evaluate(None, y0)
+            y0[: combined_submesh[0].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.zeros_like(combined_submesh[0].nodes),
         )
 
     def test_process_model_concatenation(self):
@@ -461,7 +463,7 @@ class TestDiscretise(unittest.TestCase):
 
         disc.process_model(model)
         y0 = model.concatenated_initial_conditions
-        np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh.nodes))
+        np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh[0].nodes))
 
         # grad and div are identity operators here
         np.testing.assert_array_equal(y0, model.concatenated_rhs.evaluate(None, y0))
@@ -483,7 +485,7 @@ class TestDiscretise(unittest.TestCase):
         broad = disc._spatial_methods[whole_cell[0]].broadcast(a, whole_cell)
         self.assertIsInstance(broad, pybamm.Array)
         np.testing.assert_array_equal(
-            broad.evaluate(), 7 * np.ones_like(combined_submesh.nodes)
+            broad.evaluate(), 7 * np.ones_like(combined_submesh[0].nodes)
         )
         self.assertEqual(broad.domain, whole_cell)
 
@@ -526,8 +528,8 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(eqn_disc, pybamm.Vector)
         expected_vector = np.concatenate(
             [
-                5 * np.ones_like(mesh["negative electrode"].nodes),
-                4 * np.ones_like(mesh["positive electrode"].nodes),
+                5 * np.ones_like(mesh["negative electrode"][0].nodes),
+                4 * np.ones_like(mesh["positive electrode"][0].nodes),
             ]
         )
         np.testing.assert_allclose(eqn_disc.evaluate(), expected_vector)
@@ -541,7 +543,7 @@ class TestDiscretise(unittest.TestCase):
         x1_disc = disc.process_symbol(x1)
         self.assertIsInstance(x1_disc, pybamm.Vector)
         np.testing.assert_array_equal(
-            x1_disc.evaluate(), disc.mesh["negative electrode"].nodes
+            x1_disc.evaluate(), disc.mesh["negative electrode"][0].nodes
         )
 
         x2 = pybamm.SpatialVariable("x", ["negative electrode", "separator"])
@@ -549,14 +551,14 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(x2_disc, pybamm.Vector)
         np.testing.assert_array_equal(
             x2_disc.evaluate(),
-            disc.mesh.combine_submeshes("negative electrode", "separator").nodes,
+            disc.mesh.combine_submeshes("negative electrode", "separator")[0].nodes,
         )
 
         r = 3 * pybamm.SpatialVariable("r", ["negative particle"])
         r_disc = disc.process_symbol(r)
         self.assertIsInstance(r_disc.children[1], pybamm.Vector)
         np.testing.assert_array_equal(
-            r_disc.evaluate(), 3 * disc.mesh["negative particle"].nodes
+            r_disc.evaluate(), 3 * disc.mesh["negative particle"][0].nodes
         )
 
 
diff --git a/tests/test_expression_tree/test_broadcasts.py b/tests/test_expression_tree/test_broadcasts.py
index ed4828e3..a83047b8 100644
--- a/tests/test_expression_tree/test_broadcasts.py
+++ b/tests/test_expression_tree/test_broadcasts.py
@@ -37,7 +37,7 @@ class TestBroadcasts(unittest.TestCase):
         disc = get_discretisation_for_testing()
         mesh = disc.mesh
         for dom in mesh.keys():
-            mesh[dom].npts_for_broadcast = mesh[dom].npts
+            mesh[dom][0].npts_for_broadcast = mesh[dom][0].npts
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         combined_submeshes = mesh.combine_submeshes(*whole_cell)
 
@@ -45,7 +45,7 @@ class TestBroadcasts(unittest.TestCase):
         a = pybamm.Scalar(7)
         broad = pybamm.NumpyBroadcast(a, whole_cell, mesh)
         np.testing.assert_array_equal(
-            broad.evaluate(), 7 * np.ones_like(combined_submeshes.nodes)
+            broad.evaluate(), 7 * np.ones_like(combined_submeshes[0].nodes)
         )
         self.assertEqual(broad.domain, whole_cell)
 
@@ -54,13 +54,13 @@ class TestBroadcasts(unittest.TestCase):
         broad = pybamm.NumpyBroadcast(t, whole_cell, mesh)
         self.assertEqual(broad.domain, whole_cell)
         np.testing.assert_array_equal(
-            broad.evaluate(t=3), 13 * np.ones_like(combined_submeshes.nodes)
+            broad.evaluate(t=3), 13 * np.ones_like(combined_submeshes[0].nodes)
         )
         np.testing.assert_array_equal(
             broad.evaluate(t=np.linspace(0, 1)),
             (
                 (3 * np.linspace(0, 1) + 4)[:, np.newaxis]
-                * np.ones_like(combined_submeshes.nodes)
+                * np.ones_like(combined_submeshes[0].nodes)
             ).T,
         )
 
@@ -69,7 +69,8 @@ class TestBroadcasts(unittest.TestCase):
         broad = pybamm.NumpyBroadcast(state_vec, whole_cell, mesh)
         y = np.vstack([np.linspace(0, 1), np.linspace(0, 2)])
         np.testing.assert_array_equal(
-            broad.evaluate(y=y), (y[1:2].T * np.ones_like(combined_submeshes.nodes)).T
+            broad.evaluate(y=y),
+            (y[1:2].T * np.ones_like(combined_submeshes[0].nodes)).T,
         )
 
         # state vector - bad input
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index c1625811..2bd9b17f 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -92,14 +92,14 @@ class TestConcatenations(unittest.TestCase):
         a_dom = ["negative electrode"]
         b_dom = ["positive electrode"]
         a = pybamm.NumpyBroadcast(pybamm.Scalar(2, a_dom), a_dom, mesh)
-        b = pybamm.Vector(np.ones_like(mesh[b_dom[0]].nodes), domain=b_dom)
+        b = pybamm.Vector(np.ones_like(mesh[b_dom[0]][0].nodes), domain=b_dom)
 
         # concatenate them the "wrong" way round to check they get reordered correctly
         conc = pybamm.DomainConcatenation([b, a], mesh)
         np.testing.assert_array_equal(
             conc.evaluate(),
             np.concatenate(
-                [np.full(mesh[a_dom[0]].npts, 2), np.full(mesh[b_dom[0]].npts, 1)]
+                [np.full(mesh[a_dom[0]][0].npts, 2), np.full(mesh[b_dom[0]][0].npts, 1)]
             ),
         )
 
@@ -109,7 +109,7 @@ class TestConcatenations(unittest.TestCase):
         a = pybamm.NumpyBroadcast(pybamm.Scalar(2, a_dom), a_dom, mesh)
         b = pybamm.Vector(
             np.concatenate(
-                [np.full(mesh[b_dom[0]].npts, 1), np.full(mesh[b_dom[1]].npts, 3)]
+                [np.full(mesh[b_dom[0]][0].npts, 1), np.full(mesh[b_dom[1]][0].npts, 3)]
             ),
             domain=b_dom,
         )
@@ -119,9 +119,9 @@ class TestConcatenations(unittest.TestCase):
             conc.evaluate(),
             np.concatenate(
                 [
-                    np.full(mesh[b_dom[0]].npts, 1),
-                    np.full(mesh[a_dom[0]].npts, 2),
-                    np.full(mesh[b_dom[1]].npts, 3),
+                    np.full(mesh[b_dom[0]][0].npts, 1),
+                    np.full(mesh[a_dom[0]][0].npts, 2),
+                    np.full(mesh[b_dom[1]][0].npts, 3),
                 ]
             ),
         )
@@ -166,9 +166,9 @@ class TestConcatenations(unittest.TestCase):
             processed_conc.evaluate(),
             np.concatenate(
                 [
-                    np.ones(mesh["negative electrode"].npts),
-                    2 * np.ones(mesh["separator"].npts),
-                    3 * np.ones(mesh["positive electrode"].npts),
+                    np.ones(mesh["negative electrode"][0].npts),
+                    2 * np.ones(mesh["separator"][0].npts),
+                    3 * np.ones(mesh["positive electrode"][0].npts),
                 ]
             ),
         )
@@ -191,9 +191,9 @@ class TestConcatenations(unittest.TestCase):
             processed_conc.evaluate(t=2),
             np.concatenate(
                 [
-                    2 * np.ones(mesh["negative electrode"].npts),
-                    4 * np.ones(mesh["separator"].npts),
-                    6 * np.ones(mesh["positive electrode"].npts),
+                    2 * np.ones(mesh["negative electrode"][0].npts),
+                    4 * np.ones(mesh["separator"][0].npts),
+                    6 * np.ones(mesh["positive electrode"][0].npts),
                 ]
             ),
         )
@@ -217,9 +217,9 @@ class TestConcatenations(unittest.TestCase):
             processed_conc.evaluate(y=y),
             np.concatenate(
                 [
-                    np.ones(mesh["negative electrode"].npts),
-                    2 * np.ones(mesh["separator"].npts),
-                    3 * np.ones(mesh["positive electrode"].npts),
+                    np.ones(mesh["negative electrode"][0].npts),
+                    2 * np.ones(mesh["separator"][0].npts),
+                    3 * np.ones(mesh["positive electrode"][0].npts),
                 ]
             ),
         )
@@ -239,9 +239,9 @@ class TestConcatenations(unittest.TestCase):
             processed_conc.evaluate(t=2, y=y),
             np.concatenate(
                 [
-                    np.ones(mesh["negative electrode"].npts),
-                    4 * np.ones(mesh["separator"].npts),
-                    3 * np.ones(mesh["positive electrode"].npts),
+                    np.ones(mesh["negative electrode"][0].npts),
+                    4 * np.ones(mesh["separator"][0].npts),
+                    3 * np.ones(mesh["positive electrode"][0].npts),
                 ]
             ),
         )
diff --git a/tests/test_geometry/test_geometry.py b/tests/test_geometry/test_geometry.py
index cbdcd133..3081caf2 100644
--- a/tests/test_geometry/test_geometry.py
+++ b/tests/test_geometry/test_geometry.py
@@ -22,7 +22,8 @@ class TestGeometry1DMacro(unittest.TestCase):
 
     def test_geometry_keys(self):
         geometry = pybamm.Geometry1DMacro()
-        for spatial_vars in geometry.values():
+        for prim_sec_vars in geometry.values():
+            spatial_vars = prim_sec_vars["primary"]
             all(
                 self.assertIsInstance(spatial_var, pybamm.SpatialVariable)
                 for spatial_var in spatial_vars.keys()
@@ -34,10 +35,9 @@ class TestGeometry1DMicro(unittest.TestCase):
         geometry = pybamm.Geometry1DMicro()
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         x = pybamm.SpatialVariable("x", whole_cell)
-        custom_geometry = {
-            "negative electrode": {
-                x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
-            }
+        custom_geometry = {}
+        custom_geometry["negative electrode"] = {
+            "primary": {x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}}
         }
         geometry.update(custom_geometry)
 
@@ -47,11 +47,12 @@ class TestGeometry1DMicro(unittest.TestCase):
 
     def test_geometry_keys(self):
         geometry = pybamm.Geometry1DMicro()
-        for spatial_vars in geometry.values():
-            all(
-                self.assertIsInstance(spatial_var, pybamm.SpatialVariable)
-                for spatial_var in spatial_vars.keys()
-            )
+        for prim_sec_vars in geometry.values():
+            for spatial_vars in prim_sec_vars.values():
+                all(
+                    self.assertIsInstance(spatial_var, pybamm.SpatialVariable)
+                    for spatial_var in spatial_vars.keys()
+                )
 
 
 class TestGeometry3DMacro(unittest.TestCase):
@@ -61,7 +62,7 @@ class TestGeometry3DMacro(unittest.TestCase):
         x = pybamm.SpatialVariable("x", whole_cell)
         custom_geometry = {
             "negative electrode": {
-                x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
+                "primary": {x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}}
             }
         }
 
@@ -72,11 +73,12 @@ class TestGeometry3DMacro(unittest.TestCase):
 
     def test_geometry_keys(self):
         geometry = pybamm.Geometry3DMacro()
-        for spatial_vars in geometry.values():
-            all(
-                self.assertIsInstance(spatial_var, pybamm.SpatialVariable)
-                for spatial_var in spatial_vars.keys()
-            )
+        for prim_sec_vars in geometry.values():
+            for spatial_vars in prim_sec_vars.values():
+                all(
+                    self.assertIsInstance(spatial_var, pybamm.SpatialVariable)
+                    for spatial_var in spatial_vars.keys()
+                )
 
 
 class TestGeometry(unittest.TestCase):
@@ -100,7 +102,7 @@ class TestGeometry(unittest.TestCase):
         x = pybamm.SpatialVariable("x", whole_cell)
         custom_geometry = {
             "negative electrode": {
-                x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
+                "primary": {x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}}
             }
         }
         geometry = pybamm.Geometry(
diff --git a/tests/test_meshes/test_meshes.py b/tests/test_meshes/test_meshes.py
index 9f6b2483..349a8437 100644
--- a/tests/test_meshes/test_meshes.py
+++ b/tests/test_meshes/test_meshes.py
@@ -41,18 +41,18 @@ class TestMesh(unittest.TestCase):
         mesh = mesh_type(geometry, submesh_types, submesh_pts)
 
         # check boundary locations
-        self.assertEqual(mesh["negative electrode"].edges[0], 0)
-        self.assertEqual(mesh["positive electrode"].edges[-1], 1)
+        self.assertEqual(mesh["negative electrode"][0].edges[0], 0)
+        self.assertEqual(mesh["positive electrode"][0].edges[-1], 1)
 
         # check internal boundary locations
         self.assertEqual(
-            mesh["negative electrode"].edges[-1], mesh["separator"].edges[0]
+            mesh["negative electrode"][0].edges[-1], mesh["separator"][0].edges[0]
         )
         self.assertEqual(
-            mesh["positive electrode"].edges[0], mesh["separator"].edges[-1]
+            mesh["positive electrode"][0].edges[0], mesh["separator"][0].edges[-1]
         )
         for domain in mesh:
-            self.assertEqual(len(mesh[domain].edges), len(mesh[domain].nodes) + 1)
+            self.assertEqual(len(mesh[domain][0].edges), len(mesh[domain][0].nodes) + 1)
 
     def test_mesh_sizes(self):
         param = pybamm.ParameterValues(
@@ -87,8 +87,8 @@ class TestMesh(unittest.TestCase):
         # create mesh
         mesh = mesh_type(geometry, submesh_types, submesh_pts)
         for domain in mesh:
-            self.assertEqual(mesh[domain].npts, submesh_pts[domain]["x"])
-            self.assertEqual(len(mesh[domain].edges) - 1, submesh_pts[domain]["x"])
+            self.assertEqual(mesh[domain][0].npts, submesh_pts[domain]["x"])
+            self.assertEqual(len(mesh[domain][0].edges) - 1, submesh_pts[domain]["x"])
 
     def test_combine_submeshes(self):
         param = pybamm.ParameterValues(
@@ -125,13 +125,13 @@ class TestMesh(unittest.TestCase):
 
         # create submesh
         submesh = mesh.combine_submeshes("negative electrode", "separator")
-        self.assertEqual(submesh.edges[0], 0)
-        self.assertEqual(submesh.edges[-1], mesh["separator"].edges[-1])
+        self.assertEqual(submesh[0].edges[0], 0)
+        self.assertEqual(submesh[0].edges[-1], mesh["separator"][0].edges[-1])
         self.assertAlmostEqual(
             np.linalg.norm(
-                submesh.nodes
+                submesh[0].nodes
                 - np.concatenate(
-                    [mesh["negative electrode"].nodes, mesh["separator"].nodes]
+                    [mesh["negative electrode"][0].nodes, mesh["separator"][0].nodes]
                 )
             ),
             0,
@@ -175,18 +175,54 @@ class TestMesh(unittest.TestCase):
         mesh.add_ghost_meshes()
 
         np.testing.assert_array_equal(
-            mesh["negative electrode_left ghost cell"].edges[1],
-            mesh["negative electrode"].edges[0],
+            mesh["negative electrode_left ghost cell"][0].edges[1],
+            mesh["negative electrode"][0].edges[0],
         )
         np.testing.assert_array_equal(
-            mesh["negative electrode_left ghost cell"].edges[0],
-            -mesh["negative electrode"].edges[1],
+            mesh["negative electrode_left ghost cell"][0].edges[0],
+            -mesh["negative electrode"][0].edges[1],
         )
         np.testing.assert_array_equal(
-            mesh["positive electrode_right ghost cell"].edges[0],
-            mesh["positive electrode"].edges[-1],
+            mesh["positive electrode_right ghost cell"][0].edges[0],
+            mesh["positive electrode"][0].edges[-1],
         )
 
+    def test_multiple_meshes(self):
+        param = pybamm.ParameterValues(
+            base_parameters={
+                "Negative electrode width": 0.1,
+                "Separator width": 0.2,
+                "Positive electrode width": 0.3,
+            }
+        )
+
+        geometry = pybamm.Geometry("1+1D micro")
+        param.process_geometry(geometry)
+
+        # provide mesh properties
+        submesh_pts = {
+            "negative particle": {"r": 5, "x": 10},
+            "positive particle": {"r": 6, "x": 10},
+        }
+        submesh_types = {
+            "negative particle": pybamm.Uniform1DSubMesh,
+            "positive particle": pybamm.Uniform1DSubMesh,
+        }
+
+        mesh = pybamm.Mesh(geometry, submesh_types, submesh_pts)
+
+        # check types
+        self.assertIsInstance(mesh["negative particle"], list)
+        self.assertIsInstance(mesh["positive particle"], list)
+        self.assertEqual(len(mesh["negative particle"]), 10)
+        self.assertEqual(len(mesh["positive particle"]), 10)
+
+        for i in range(10):
+            self.assertIsInstance(mesh["negative particle"][i], pybamm.Uniform1DSubMesh)
+            self.assertIsInstance(mesh["positive particle"][i], pybamm.Uniform1DSubMesh)
+            self.assertEqual(mesh["negative particle"][i].npts, 5)
+            self.assertEqual(mesh["positive particle"][i].npts, 6)
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_models/test_interface.py b/tests/test_models/test_interface.py
index 5e6b053b..02d0a7ae 100644
--- a/tests/test_models/test_interface.py
+++ b/tests/test_models/test_interface.py
@@ -80,15 +80,14 @@ class TestHomogeneousReaction(unittest.TestCase):
 
         submesh = disc.mesh.combine_submeshes(*whole_cell)
 
-        # processed_rxn should be a concatenation with the right shape
         self.assertIsInstance(processed_rxn, pybamm.Concatenation)
-        self.assertEqual(processed_rxn.evaluate(0, None).shape, submesh.nodes.shape)
+        self.assertEqual(processed_rxn.evaluate(0, None).shape, submesh[0].nodes.shape)
 
         # test values
         l_n = param.process_symbol(pybamm.standard_parameters.l_n)
         l_p = param.process_symbol(pybamm.standard_parameters.l_p)
-        npts_n = mesh["negative electrode"].npts
-        npts_s = mesh["separator"].npts
+        npts_n = mesh["negative electrode"][0].npts
+        npts_s = mesh["separator"][0].npts
         np.testing.assert_array_equal(
             (l_n * processed_rxn).evaluate(0, None)[:npts_n], 1
         )
@@ -315,18 +314,18 @@ class TestButlerVolmer(unittest.TestCase):
         processed_bv_p = disc.process_symbol(param_bv_p)
 
         submesh = np.concatenate(
-            [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
+            [mesh["negative electrode"][0].nodes, mesh["positive electrode"][0].nodes]
         )
         y = np.concatenate([submesh ** 2, submesh ** 3])
 
         # should evaluate to vectors with the right shape
         self.assertEqual(
             processed_bv_n.evaluate(None, y).shape,
-            mesh["negative electrode"].nodes.shape,
+            mesh["negative electrode"][0].nodes.shape,
         )
         self.assertEqual(
             processed_bv_p.evaluate(None, y).shape,
-            mesh["positive electrode"].nodes.shape,
+            mesh["positive electrode"][0].nodes.shape,
         )
 
     def test_discretisation_with_particles(self):
@@ -370,18 +369,18 @@ class TestButlerVolmer(unittest.TestCase):
         processed_bv_p = disc.process_symbol(param_bv_p)
 
         submesh = np.concatenate(
-            [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
+            [mesh["negative electrode"][0].nodes, mesh["positive electrode"][0].nodes]
         )
         y = np.concatenate([submesh ** 2, submesh ** 3, submesh])
 
         # should evaluate to vectors with the right shape
         self.assertEqual(
             processed_bv_n.evaluate(None, y).shape,
-            mesh["negative electrode"].nodes.shape,
+            mesh["negative electrode"][0].nodes.shape,
         )
         self.assertEqual(
             processed_bv_p.evaluate(None, y).shape,
-            mesh["positive electrode"].nodes.shape,
+            mesh["positive electrode"][0].nodes.shape,
         )
 
     def test_discretisation_whole(self):
@@ -413,9 +412,11 @@ class TestButlerVolmer(unittest.TestCase):
         # test
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         whole_cell_mesh = disc.mesh.combine_submeshes(*whole_cell)
-        y = np.concatenate([whole_cell_mesh.nodes ** 2, whole_cell_mesh.nodes ** 3])
+        y = np.concatenate(
+            [whole_cell_mesh[0].nodes ** 2, whole_cell_mesh[0].nodes ** 3]
+        )
         self.assertEqual(
-            processed_bv_whole.evaluate(None, y).shape, whole_cell_mesh.nodes.shape
+            processed_bv_whole.evaluate(None, y).shape, whole_cell_mesh[0].nodes.shape
         )
 
     def test_discretisation_whole_with_particles(self):
@@ -451,13 +452,13 @@ class TestButlerVolmer(unittest.TestCase):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         whole_cell_mesh = disc.mesh.combine_submeshes(*whole_cell)
         submesh = np.concatenate(
-            [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
+            [mesh["negative electrode"][0].nodes, mesh["positive electrode"][0].nodes]
         )
         y = np.concatenate(
-            [whole_cell_mesh.nodes ** 2, whole_cell_mesh.nodes ** 3, submesh]
+            [whole_cell_mesh[0].nodes ** 2, whole_cell_mesh[0].nodes ** 3, submesh]
         )
         self.assertEqual(
-            processed_bv_whole.evaluate(None, y).shape, whole_cell_mesh.nodes.shape
+            processed_bv_whole.evaluate(None, y).shape, whole_cell_mesh[0].nodes.shape
         )
 
 
@@ -548,15 +549,15 @@ class TestExchangeCurrentDensity(unittest.TestCase):
         processed_j0p = disc.process_symbol(param_j0p)
 
         submesh = np.concatenate(
-            [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
+            [mesh["negative electrode"][0].nodes, mesh["positive electrode"][0].nodes]
         )
         y = submesh ** 2
         # should evaluate to vectors with the right shape
         self.assertEqual(
-            processed_j0n.evaluate(y=y).shape, mesh["negative electrode"].nodes.shape
+            processed_j0n.evaluate(y=y).shape, mesh["negative electrode"][0].nodes.shape
         )
         self.assertEqual(
-            processed_j0p.evaluate(y=y).shape, mesh["positive electrode"].nodes.shape
+            processed_j0p.evaluate(y=y).shape, mesh["positive electrode"][0].nodes.shape
         )
 
     def test_discretisation_surface_conc(self):
@@ -582,19 +583,19 @@ class TestExchangeCurrentDensity(unittest.TestCase):
 
         submesh = np.concatenate(
             [
-                mesh["negative electrode"].nodes,
-                mesh["positive electrode"].nodes,
-                mesh["negative electrode"].nodes,
-                mesh["positive electrode"].nodes,
+                mesh["negative electrode"][0].nodes,
+                mesh["positive electrode"][0].nodes,
+                mesh["negative electrode"][0].nodes,
+                mesh["positive electrode"][0].nodes,
             ]
         )
         y = submesh ** 2
         # should evaluate to vectors with the right shape
         self.assertEqual(
-            processed_j0n.evaluate(y=y).shape, mesh["negative electrode"].nodes.shape
+            processed_j0n.evaluate(y=y).shape, mesh["negative electrode"][0].nodes.shape
         )
         self.assertEqual(
-            processed_j0p.evaluate(y=y).shape, mesh["positive electrode"].nodes.shape
+            processed_j0p.evaluate(y=y).shape, mesh["positive electrode"][0].nodes.shape
         )
 
 
diff --git a/tests/test_models/test_simple_ode_model.py b/tests/test_models/test_simple_ode_model.py
index 20d41292..4d169484 100644
--- a/tests/test_models/test_simple_ode_model.py
+++ b/tests/test_models/test_simple_ode_model.py
@@ -32,13 +32,13 @@ class TestSimpleODEModel(unittest.TestCase):
         )
         np.testing.assert_array_almost_equal(
             model.variables["b broadcasted"].evaluate(T, Y),
-            np.ones((combined_submesh.npts, T.size)),
+            np.ones((combined_submesh[0].npts, T.size)),
         )
         np.testing.assert_array_almost_equal(
             model.variables["c broadcasted"].evaluate(T, Y),
-            np.ones(sum([mesh[d].npts for d in ["negative electrode", "separator"]]))[
-                :, np.newaxis
-            ]
+            np.ones(
+                sum([mesh[d][0].npts for d in ["negative electrode", "separator"]])
+            )[:, np.newaxis]
             * np.exp(-T),
         )
 
diff --git a/tests/test_parameters/test_standard_parameters_lead_acid.py b/tests/test_parameters/test_standard_parameters_lead_acid.py
index 1ccd4513..9afaaa15 100644
--- a/tests/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/test_parameters/test_standard_parameters_lead_acid.py
@@ -74,7 +74,7 @@ class TestStandardParametersLeadAcid(unittest.TestCase):
         combined_submeshes = disc.mesh.combine_submeshes(
             "negative electrode", "separator", "positive electrode"
         )
-        self.assertEqual(processed_s.shape, combined_submeshes.nodes.shape)
+        self.assertEqual(processed_s.shape, combined_submeshes[0].nodes.shape)
 
     def test_current_functions(self):
         # create current functions
diff --git a/tests/test_spatial_methods/test_finite_volume.py b/tests/test_spatial_methods/test_finite_volume.py
index c8c24e09..42e92b14 100644
--- a/tests/test_spatial_methods/test_finite_volume.py
+++ b/tests/test_spatial_methods/test_finite_volume.py
@@ -2,7 +2,7 @@
 # Test for the operator class
 #
 import pybamm
-from tests import get_mesh_for_testing
+from tests import get_mesh_for_testing, get_p2d_mesh_for_testing
 
 import numpy as np
 import unittest
@@ -47,12 +47,12 @@ class TestFiniteVolume(unittest.TestCase):
         surf_eqn_disc = disc.process_symbol(surf_eqn)
 
         # check constant extrapolates to constant
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         self.assertEqual(surf_eqn_disc.evaluate(None, constant_y), 1)
 
         # check linear variable extrapolates correctly
-        linear_y = combined_submesh.nodes
-        y_surf = combined_submesh.nodes[-1] + combined_submesh.d_nodes[-1] / 2
+        linear_y = combined_submesh[0].nodes
+        y_surf = combined_submesh[0].nodes[-1] + combined_submesh[0].d_nodes[-1] / 2
         self.assertEqual(surf_eqn_disc.evaluate(None, linear_y), y_surf)
 
     def test_discretise_diffusivity_times_spatial_operator(self):
@@ -69,7 +69,7 @@ class TestFiniteVolume(unittest.TestCase):
         # 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)
+        y_test = np.ones_like(combined_submesh[0].nodes)
         for eqn in [
             var * pybamm.grad(var),
             var ** 2 * pybamm.grad(var),
@@ -125,12 +125,12 @@ class TestFiniteVolume(unittest.TestCase):
         lbc = pybamm.Scalar(0)
         rbc = pybamm.Scalar(3)
         symbol_plus_ghost = pybamm.FiniteVolume(mesh).add_ghost_nodes(
-            discretised_symbol, lbc, rbc
+            var, discretised_symbol, lbc, rbc
         )
 
         # Test
         combined_submesh = mesh.combine_submeshes(*whole_cell)
-        y_test = np.ones_like(combined_submesh.nodes)
+        y_test = np.ones_like(combined_submesh[0].nodes)
         np.testing.assert_array_equal(
             symbol_plus_ghost.evaluate(None, y_test)[1:-1],
             discretised_symbol.evaluate(None, y_test),
@@ -152,6 +152,77 @@ class TestFiniteVolume(unittest.TestCase):
             3,
         )
 
+    def test_p2d_add_ghost_nodes(self):
+        # create discretisation
+        mesh = get_p2d_mesh_for_testing()
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+            "positive particle": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        # add ghost nodes
+        c_s_n = pybamm.Variable("c_s_n", domain=["negative particle"])
+        c_s_p = pybamm.Variable("c_s_p", domain=["positive particle"])
+
+        disc.set_variable_slices([c_s_n])
+        disc_c_s_n = pybamm.StateVector(disc._y_slices[c_s_n.id])
+
+        disc.set_variable_slices([c_s_p])
+        disc_c_s_p = pybamm.StateVector(disc._y_slices[c_s_p.id])
+        lbc = pybamm.Scalar(0)
+        rbc = pybamm.Scalar(3)
+        c_s_n_plus_ghost = pybamm.FiniteVolume(mesh).add_ghost_nodes(
+            c_s_n, disc_c_s_n, lbc, rbc
+        )
+        c_s_p_plus_ghost = pybamm.FiniteVolume(mesh).add_ghost_nodes(
+            c_s_p, disc_c_s_p, lbc, rbc
+        )
+
+        mesh_s_n = mesh["negative particle"]
+        mesh_s_p = mesh["positive particle"]
+
+        n_prim_pts = mesh_s_n[0].npts
+        n_sec_pts = len(mesh_s_n)
+
+        p_prim_pts = mesh_s_p[0].npts
+        p_sec_pts = len(mesh_s_p)
+
+        y_s_n_test = np.kron(np.ones(n_sec_pts), np.ones(n_prim_pts))
+        y_s_p_test = np.kron(np.ones(p_sec_pts), np.ones(p_prim_pts))
+
+        # evaluate with and without ghost points
+        c_s_n_eval = disc_c_s_n.evaluate(None, y_s_n_test)
+        c_s_n_ghost_eval = c_s_n_plus_ghost.evaluate(None, y_s_n_test)
+
+        c_s_p_eval = disc_c_s_p.evaluate(None, y_s_p_test)
+        c_s_p_ghost_eval = c_s_p_plus_ghost.evaluate(None, y_s_p_test)
+
+        # reshape to make easy to deal with
+        c_s_n_eval = np.reshape(c_s_n_eval, [n_sec_pts, n_prim_pts])
+        c_s_n_ghost_eval = np.reshape(c_s_n_ghost_eval, [n_sec_pts, n_prim_pts + 2])
+
+        c_s_p_eval = np.reshape(c_s_p_eval, [p_sec_pts, p_prim_pts])
+        c_s_p_ghost_eval = np.reshape(c_s_p_ghost_eval, [p_sec_pts, p_prim_pts + 2])
+
+        np.testing.assert_array_equal(c_s_n_ghost_eval[:, 1:-1], c_s_n_eval)
+        np.testing.assert_array_equal(c_s_p_ghost_eval[:, 1:-1], c_s_p_eval)
+
+        np.testing.assert_array_equal(
+            (c_s_n_ghost_eval[:, 0] + c_s_n_ghost_eval[:, 1]) / 2, 0
+        )
+        np.testing.assert_array_equal(
+            (c_s_p_ghost_eval[:, 0] + c_s_p_ghost_eval[:, 1]) / 2, 0
+        )
+
+        np.testing.assert_array_equal(
+            (c_s_n_ghost_eval[:, -2] + c_s_n_ghost_eval[:, -1]) / 2, 3
+        )
+        np.testing.assert_array_equal(
+            (c_s_p_ghost_eval[:, -2] + c_s_p_ghost_eval[:, -1]) / 2, 3
+        )
+
     def test_grad_div_shapes_Dirichlet_bcs(self):
         """
         Test grad and div with Dirichlet boundary conditions (applied by grad on var)
@@ -178,14 +249,14 @@ class TestFiniteVolume(unittest.TestCase):
         disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges),
+            np.zeros_like(combined_submesh[0].edges),
         )
 
         # div: test on linear y (should have laplacian zero) so change bcs
-        linear_y = combined_submesh.nodes
+        linear_y = combined_submesh[0].nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -196,12 +267,14 @@ class TestFiniteVolume(unittest.TestCase):
 
         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)
+            grad_eqn_disc.evaluate(None, linear_y),
+            np.ones_like(combined_submesh[0].edges),
         )
 
         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)
+            div_eqn_disc.evaluate(None, linear_y),
+            np.zeros_like(combined_submesh[0].nodes),
         )
 
     def test_spherical_grad_div_shapes_Dirichlet_bcs(self):
@@ -231,10 +304,10 @@ class TestFiniteVolume(unittest.TestCase):
         disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges),
+            np.zeros_like(combined_submesh[0].edges),
         )
 
         boundary_conditions = {
@@ -242,15 +315,16 @@ class TestFiniteVolume(unittest.TestCase):
         }
         disc._bcs = boundary_conditions
 
-        y_linear = combined_submesh.nodes
+        y_linear = combined_submesh[0].nodes
         grad_eqn_disc = disc.process_symbol(grad_eqn)
         np.testing.assert_array_almost_equal(
-            grad_eqn_disc.evaluate(None, y_linear), np.ones_like(combined_submesh.edges)
+            grad_eqn_disc.evaluate(None, y_linear),
+            np.ones_like(combined_submesh[0].edges),
         )
 
         # div: test on linear r^2
         # div (grad r^2) = 6
-        const = 6 * np.ones(combined_submesh.npts)
+        const = 6 * np.ones(combined_submesh[0].npts)
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
@@ -260,9 +334,63 @@ class TestFiniteVolume(unittest.TestCase):
 
         div_eqn_disc = disc.process_symbol(div_eqn)
         np.testing.assert_array_almost_equal(
-            div_eqn_disc.evaluate(None, const), np.zeros_like(combined_submesh.nodes)
+            div_eqn_disc.evaluate(None, const), np.zeros_like(combined_submesh[0].nodes)
         )
 
+    def test_p2d_spherical_grad_div_shapes_Dirichlet_bcs(self):
+        """
+        Test grad and div with Dirichlet boundary conditions (applied by grad on var)
+        in the pseudo 2-dimensional case
+        """
+
+        mesh = get_p2d_mesh_for_testing()
+        spatial_methods = {
+            "macroscale": pybamm.FiniteVolume,
+            "negative particle": pybamm.FiniteVolume,
+            "positive particle": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        n_mesh = mesh["negative particle"]
+
+        mesh.add_ghost_meshes()
+        disc.mesh.add_ghost_meshes()
+
+        var = pybamm.Variable("var", domain=["negative particle"])
+        grad_eqn = pybamm.grad(var)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
+        }
+        disc._bcs = boundary_conditions
+
+        disc.set_variable_slices([var])
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+
+        prim_pts = n_mesh[0].npts
+        sec_pts = len(n_mesh)
+        constant_y = np.kron(np.ones(sec_pts), np.ones(prim_pts))
+
+        grad_eval = grad_eqn_disc.evaluate(None, constant_y)
+        grad_eval = np.reshape(grad_eval, [sec_pts, prim_pts + 1])
+
+        np.testing.assert_array_equal(grad_eval, np.zeros([sec_pts, prim_pts + 1]))
+
+        # div: test on linear r^2
+        # div (grad r^2) = 6
+        const = 6 * np.ones(sec_pts * prim_pts)
+
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(6), "right": pybamm.Scalar(6)}
+        }
+        disc._bcs = boundary_conditions
+
+        div_eqn_disc = disc.process_symbol(div_eqn)
+        div_eval = div_eqn_disc.evaluate(None, const)
+        div_eval = np.reshape(div_eval, [sec_pts, prim_pts])
+        np.testing.assert_array_almost_equal(div_eval, np.zeros([sec_pts, prim_pts]))
+
     def test_grad_div_shapes_Neumann_bcs(self):
         """Test grad and div with Neumann boundary conditions (applied by div on N)"""
         whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -280,10 +408,10 @@ class TestFiniteVolume(unittest.TestCase):
         disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[1:-1]),
+            np.zeros_like(combined_submesh[0].edges[1:-1]),
         )
 
         # div
@@ -296,13 +424,14 @@ class TestFiniteVolume(unittest.TestCase):
         div_eqn_disc = disc.process_symbol(div_eqn)
 
         # Linear y should have laplacian zero
-        linear_y = combined_submesh.nodes
+        linear_y = combined_submesh[0].nodes
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[1:-1]),
+            np.ones_like(combined_submesh[0].edges[1:-1]),
         )
         np.testing.assert_array_almost_equal(
-            div_eqn_disc.evaluate(None, linear_y), np.zeros_like(combined_submesh.nodes)
+            div_eqn_disc.evaluate(None, linear_y),
+            np.zeros_like(combined_submesh[0].nodes),
         )
 
     def test_spherical_grad_div_shapes_Neumann_bcs(self):
@@ -321,16 +450,16 @@ class TestFiniteVolume(unittest.TestCase):
         disc.set_variable_slices([var])
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[1:-1]),
+            np.zeros_like(combined_submesh[0].edges[1:-1]),
         )
 
-        linear_y = combined_submesh.nodes
+        linear_y = combined_submesh[0].nodes
         np.testing.assert_array_almost_equal(
             grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[1:-1]),
+            np.ones_like(combined_submesh[0].edges[1:-1]),
         )
         # div
         # div ( grad(r^2) ) == 6 , N_left = N_right = 0
@@ -342,13 +471,58 @@ class TestFiniteVolume(unittest.TestCase):
         disc._bcs = boundary_conditions
         div_eqn_disc = disc.process_symbol(div_eqn)
 
-        linear_y = combined_submesh.nodes
-        const = 6 * np.ones(combined_submesh.npts)
+        linear_y = combined_submesh[0].nodes
+        const = 6 * np.ones(combined_submesh[0].npts)
 
         np.testing.assert_array_almost_equal(
-            div_eqn_disc.evaluate(None, const), np.zeros_like(combined_submesh.nodes)
+            div_eqn_disc.evaluate(None, const), np.zeros_like(combined_submesh[0].nodes)
         )
 
+    def test_p2d_spherical_grad_div_shapes_Neumann_bcs(self):
+        """
+        Test grad and div with Dirichlet boundary conditions (applied by grad on var)
+        in the pseudo 2-dimensional case
+        """
+
+        mesh = get_p2d_mesh_for_testing()
+        spatial_methods = {"negative particle": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(mesh, spatial_methods)
+
+        n_mesh = mesh["negative particle"]
+
+        mesh.add_ghost_meshes()
+        disc.mesh.add_ghost_meshes()
+
+        # test grad
+        var = pybamm.Variable("var", domain=["negative particle"])
+        grad_eqn = pybamm.grad(var)
+        disc.set_variable_slices([var])
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
+
+        prim_pts = n_mesh[0].npts
+        sec_pts = len(n_mesh)
+        constant_y = np.kron(np.ones(sec_pts), np.ones(prim_pts))
+
+        grad_eval = grad_eqn_disc.evaluate(None, constant_y)
+        grad_eval = np.reshape(grad_eval, [sec_pts, prim_pts - 1])
+
+        np.testing.assert_array_equal(grad_eval, np.zeros([sec_pts, prim_pts - 1]))
+
+        # div
+        # div (grad r^2) = 6, N_left = N_right = 0
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            N.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
+        }
+        disc._bcs = boundary_conditions
+        div_eqn_disc = disc.process_symbol(div_eqn)
+
+        const = 6 * np.ones(sec_pts * prim_pts)
+        div_eval = div_eqn_disc.evaluate(None, const)
+        div_eval = np.reshape(div_eval, [sec_pts, prim_pts])
+        np.testing.assert_array_almost_equal(div_eval, np.zeros([sec_pts, prim_pts]))
+
     def test_grad_div_shapes_mixed_domain(self):
         """
         Test grad and div with Dirichlet boundary conditions (applied by grad on var)
@@ -374,32 +548,34 @@ class TestFiniteVolume(unittest.TestCase):
         grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         np.testing.assert_array_equal(
             grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges),
+            np.zeros_like(combined_submesh[0].edges),
         )
 
         # div: test on linear y (should have laplacian zero) so change bcs
-        linear_y = combined_submesh.nodes
+        linear_y = combined_submesh[0].nodes
         N = pybamm.grad(var)
         div_eqn = pybamm.div(N)
         boundary_conditions = {
             var.id: {
                 "left": pybamm.Scalar(0),
-                "right": pybamm.Scalar(combined_submesh.edges[-1]),
+                "right": pybamm.Scalar(combined_submesh[0].edges[-1]),
             }
         }
         disc._bcs = boundary_conditions
 
         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)
+            grad_eqn_disc.evaluate(None, linear_y),
+            np.ones_like(combined_submesh[0].edges),
         )
 
         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)
+            div_eqn_disc.evaluate(None, linear_y),
+            np.zeros_like(combined_submesh[0].nodes),
         )
 
     def test_definite_integral(self):
@@ -412,8 +588,8 @@ class TestFiniteVolume(unittest.TestCase):
         }
         disc = pybamm.Discretisation(mesh, spatial_methods)
         # lengths
-        ln = mesh["negative electrode"].edges[-1]
-        ls = mesh["separator"].edges[-1] - ln
+        ln = mesh["negative electrode"][0].edges[-1]
+        ls = mesh["separator"][0].edges[-1] - ln
         lp = 1 - (ln + ls)
 
         # macroscale variable
@@ -424,13 +600,13 @@ class TestFiniteVolume(unittest.TestCase):
         integral_eqn_disc = disc.process_symbol(integral_eqn)
 
         combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), ln + ls)
-        linear_y = combined_submesh.nodes
+        linear_y = combined_submesh[0].nodes
         self.assertAlmostEqual(
             integral_eqn_disc.evaluate(None, linear_y), (ln + ls) ** 2 / 2
         )
-        cos_y = np.cos(combined_submesh.nodes)
+        cos_y = np.cos(combined_submesh[0].nodes)
         self.assertAlmostEqual(
             integral_eqn_disc.evaluate(None, cos_y), np.sin(ln + ls), places=4
         )
@@ -443,13 +619,13 @@ class TestFiniteVolume(unittest.TestCase):
         integral_eqn_disc = disc.process_symbol(integral_eqn)
 
         combined_submesh = mesh.combine_submeshes("separator", "positive electrode")
-        constant_y = np.ones_like(combined_submesh.nodes)
+        constant_y = np.ones_like(combined_submesh[0].nodes)
         self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), ls + lp)
-        linear_y = combined_submesh.nodes
+        linear_y = combined_submesh[0].nodes
         self.assertAlmostEqual(
             integral_eqn_disc.evaluate(None, linear_y), (1 - (ln) ** 2) / 2
         )
-        cos_y = np.cos(combined_submesh.nodes)
+        cos_y = np.cos(combined_submesh[0].nodes)
         self.assertAlmostEqual(
             integral_eqn_disc.evaluate(None, cos_y), np.sin(1) - np.sin(ln), places=4
         )
@@ -461,13 +637,13 @@ class TestFiniteVolume(unittest.TestCase):
         disc.set_variable_slices([var])
         integral_eqn_disc = disc.process_symbol(integral_eqn)
 
-        constant_y = np.ones_like(mesh["negative particle"].nodes)
+        constant_y = np.ones_like(mesh["negative particle"][0].nodes)
         self.assertEqual(integral_eqn_disc.evaluate(None, constant_y), np.pi)
-        linear_y = mesh["negative particle"].nodes
+        linear_y = mesh["negative particle"][0].nodes
         self.assertAlmostEqual(
             integral_eqn_disc.evaluate(None, linear_y), 2 * np.pi / 3, places=4
         )
-        one_over_y = 1 / mesh["negative particle"].nodes
+        one_over_y = 1 / mesh["negative particle"][0].nodes
         self.assertEqual(integral_eqn_disc.evaluate(None, one_over_y), 2 * np.pi)
 
     @unittest.skip("indefinite integral not yet implemented")
@@ -565,8 +741,8 @@ class TestFiniteVolume(unittest.TestCase):
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
             # Define exact solutions
-            y = np.sin(combined_submesh.nodes)
-            grad_exact = np.cos(combined_submesh.edges[1:-1])
+            y = np.sin(combined_submesh[0].nodes)
+            grad_exact = np.cos(combined_submesh[0].edges[1:-1])
 
             # Discretise and evaluate
             disc.set_variable_slices([var])
@@ -612,8 +788,8 @@ class TestFiniteVolume(unittest.TestCase):
             disc.mesh.add_ghost_meshes()
 
             # Define exact solutions
-            y = np.sin(combined_submesh.nodes)
-            grad_exact = np.cos(combined_submesh.edges)
+            y = np.sin(combined_submesh[0].nodes)
+            grad_exact = np.cos(combined_submesh[0].edges)
 
             # Discretise and evaluate
             disc.set_variable_slices([var])
@@ -655,8 +831,8 @@ class TestFiniteVolume(unittest.TestCase):
             combined_submesh = mesh.combine_submeshes(*whole_cell)
 
             # Define exact solutions
-            y = np.sin(combined_submesh.nodes)
-            div_exact_internal = -np.sin(combined_submesh.nodes[1:-1])
+            y = np.sin(combined_submesh[0].nodes)
+            div_exact_internal = -np.sin(combined_submesh[0].nodes[1:-1])
 
             # Discretise and evaluate
             disc.set_variable_slices([var])
@@ -699,8 +875,8 @@ class TestFiniteVolume(unittest.TestCase):
             combined_submesh = mesh.combine_submeshes(*whole_cell)
 
             # Define exact solutions
-            y = np.sin(combined_submesh.nodes)
-            div_exact_internal = -np.sin(combined_submesh.nodes[1:-1])
+            y = np.sin(combined_submesh[0].nodes)
+            div_exact_internal = -np.sin(combined_submesh[0].nodes[1:-1])
 
             # Discretise and evaluate
             disc.set_variable_slices([var])
@@ -738,7 +914,7 @@ class TestFiniteVolume(unittest.TestCase):
             disc = pybamm.Discretisation(mesh, spatial_methods)
             disc._bcs = boundary_conditions
             mesh = disc.mesh["negative particle"]
-            r = mesh.nodes
+            r = mesh[0].nodes
 
             # exact solution
             y = np.sin(r)
@@ -765,6 +941,56 @@ class TestFiniteVolume(unittest.TestCase):
         rates = np.log2(errs[:-1] / errs[1:])
         np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
 
+    def test_p2d_spherical_operators(self):
+        # test div( grad( sin(r) )) == (2/r)*cos(r) - *sin(r)
+
+        domain = ["negative particle"]
+        c = pybamm.Variable("c", domain=domain)
+        N = pybamm.grad(c)
+        eqn = pybamm.div(N)
+        boundary_conditions = {
+            N.id: {"left": pybamm.Scalar(np.cos(0)), "right": pybamm.Scalar(np.cos(1))}
+        }
+
+        def get_l2_error(m):
+            mesh = get_p2d_mesh_for_testing(3, m)
+            spatial_methods = {"negative particle": pybamm.FiniteVolume}
+            disc = pybamm.Discretisation(mesh, spatial_methods)
+            disc._bcs = boundary_conditions
+            mesh = disc.mesh["negative particle"]
+            r = mesh[0].nodes
+
+            prim_pts = mesh[0].npts
+            sec_pts = len(mesh)
+
+            # exact solution
+            y = np.kron(np.ones(sec_pts), np.sin(r))
+            exact = (2 / r) * np.cos(r) - np.sin(r)
+            exact_internal = np.kron(np.ones(sec_pts), exact[1:-1])
+
+            # discretise and evaluate
+            variables = [c]
+            disc.set_variable_slices(variables)
+            eqn_disc = disc.process_symbol(eqn)
+            approx_eval = eqn_disc.evaluate(None, y)
+            approx_eval = np.reshape(approx_eval, [sec_pts, prim_pts])
+            approx_internal = approx_eval[:, 1:-1]
+            approx_internal = np.reshape(approx_internal, [sec_pts * (prim_pts - 2)])
+
+            # error
+            error = np.linalg.norm(approx_internal - exact_internal) / np.linalg.norm(
+                exact_internal
+            )
+            return error
+
+        # Get errors
+        ns = 10 * (2 ** np.arange(2, 7))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
+
+        # Get rates: expect h**1.5 convergence because of boundary conditions
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
+
     def test_discretise_spatial_variable(self):
         # create discretisation
         mesh = get_mesh_for_testing()
@@ -780,7 +1006,7 @@ class TestFiniteVolume(unittest.TestCase):
         x1_disc = disc.process_symbol(x1)
         self.assertIsInstance(x1_disc, pybamm.Vector)
         np.testing.assert_array_equal(
-            x1_disc.evaluate(), disc.mesh["negative electrode"].nodes
+            x1_disc.evaluate(), disc.mesh["negative electrode"][0].nodes
         )
 
         z = pybamm.SpatialVariable("z", ["negative electrode"])
@@ -792,14 +1018,14 @@ class TestFiniteVolume(unittest.TestCase):
         self.assertIsInstance(x2_disc, pybamm.Vector)
         np.testing.assert_array_equal(
             x2_disc.evaluate(),
-            disc.mesh.combine_submeshes("negative electrode", "separator").nodes,
+            disc.mesh.combine_submeshes("negative electrode", "separator")[0].nodes,
         )
 
         r = 3 * pybamm.SpatialVariable("r", ["negative particle"])
         r_disc = disc.process_symbol(r)
         self.assertIsInstance(r_disc.children[1], pybamm.Vector)
         np.testing.assert_array_equal(
-            r_disc.evaluate(), 3 * disc.mesh["negative particle"].nodes
+            r_disc.evaluate(), 3 * disc.mesh["negative particle"][0].nodes
         )