diff --git a/.flake8 b/.flake8
index 70990935..9acc35b7 100644
--- a/.flake8
+++ b/.flake8
@@ -5,6 +5,12 @@ exclude=
     problems,
     __init__.py,
     venv,
+    bin,
+    etc,
+    lib,
+    lib64,
+    share,
+    pyvenv.cfg
 ignore=
     # False positive for white space before ':' on list slice
     # black should format these correctly
diff --git a/.gitignore b/.gitignore
index 65bb1425..4db0408e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -46,3 +46,9 @@ pyproject.toml
 
 # virtual enviroment
 venv/
+bin/
+etc/
+lib/
+lib64
+share/
+pyvenv.cfg
diff --git a/docs/source/discretisations/base_discretisation.rst b/docs/source/discretisations/base_discretisation.rst
new file mode 100644
index 00000000..95531920
--- /dev/null
+++ b/docs/source/discretisations/base_discretisation.rst
@@ -0,0 +1,5 @@
+Base Discretisation
+===================
+
+.. autoclass:: pybamm.BaseDiscretisation
+  :members:
diff --git a/docs/source/discretisations/core.rst b/docs/source/discretisations/core.rst
deleted file mode 100644
index 1306f77e..00000000
--- a/docs/source/discretisations/core.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-Core Discretisation
-===================
-
-.. autoclass:: pybamm.BaseDiscretisation
-  :members:
-
-.. autoclass:: pybamm.MatrixVectorDiscretisation
-  :members:
diff --git a/docs/source/discretisations/index.rst b/docs/source/discretisations/index.rst
index 50c9e2bc..d25791de 100644
--- a/docs/source/discretisations/index.rst
+++ b/docs/source/discretisations/index.rst
@@ -3,6 +3,6 @@ Discretisations
 
 .. toctree::
 
-  core
+  base_discretisation
   finite_volume_discretisations
   finite_volume_meshes
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index c392e99c..59856547 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -102,10 +102,7 @@ from .parameters import standard_parameters
 #
 # Mesh and Discretisation classes
 #
-from .discretisations.base_discretisation import (
-    BaseDiscretisation,
-    MatrixVectorDiscretisation,
-)
+from .discretisations.base_discretisation import BaseDiscretisation
 from .discretisations.finite_volume_discretisations import FiniteVolumeDiscretisation
 from .discretisations.base_mesh import KNOWN_DOMAINS
 from .discretisations.base_mesh import BaseMesh, BaseSubmesh
diff --git a/pybamm/discretisations/base_discretisation.py b/pybamm/discretisations/base_discretisation.py
index 7ae9efde..a19743e4 100644
--- a/pybamm/discretisations/base_discretisation.py
+++ b/pybamm/discretisations/base_discretisation.py
@@ -294,78 +294,3 @@ class BaseDiscretisation(object):
 
     def concatenate(self, *symbols):
         return pybamm.NumpyConcatenation(*symbols)
-
-
-class MatrixVectorDiscretisation(BaseDiscretisation):
-    """
-    A base class for discretisations where the gradient and divergence operators are
-    always matrix-vector multiplications.
-
-    Parameters
-    ----------
-     mesh : :class:`pybamm.BaseMesh` (or subclass)
-            The underlying mesh for discretisation
-
-    **Extends:** :class:`pybamm.BaseDiscretisation`
-    """
-
-    def __init__(self, mesh):
-        super().__init__(mesh)
-
-    def gradient(self, symbol, y_slices, boundary_conditions):
-        """Matrix-vector multiplication to implement the gradient operator.
-        See :meth:`pybamm.BaseDiscretisation.gradient`
-        """
-        # Check that boundary condition keys are hashes (ids)
-        for key in boundary_conditions.keys():
-            assert isinstance(key, int), TypeError(
-                "boundary condition keys should be hashes, not {}".format(type(key))
-            )
-        # Discretise symbol
-        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
-        # Add boundary conditions if defined
-        if symbol.id in boundary_conditions:
-            lbc = boundary_conditions[symbol.id]["left"]
-            rbc = boundary_conditions[symbol.id]["right"]
-            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
-        gradient_matrix = self.gradient_matrix(symbol.domain)
-        return gradient_matrix * discretised_symbol
-
-    def gradient_matrix(self, domain):
-        """The gradient matrix
-
-        Parameters
-        ----------
-        domain : list
-            The domain(s) in which to compute the gradient matrix
-        """
-        raise NotImplementedError
-
-    def divergence(self, symbol, y_slices, boundary_conditions):
-        """Matrix-vector multiplication to implement the divergence operator.
-        See :meth:`pybamm.BaseDiscretisation.gradient`
-        """
-        # Check that boundary condition keys are hashes (ids)
-        for key in boundary_conditions.keys():
-            assert isinstance(key, int), TypeError(
-                "boundary condition keys should be hashes, not {}".format(type(key))
-            )
-        # Discretise symbol
-        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
-        # Add boundary conditions if defined
-        if symbol.id in boundary_conditions:
-            lbc = boundary_conditions[symbol.id]["left"]
-            rbc = boundary_conditions[symbol.id]["right"]
-            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
-        divergence_matrix = self.divergence_matrix(symbol.domain)
-        return divergence_matrix * discretised_symbol
-
-    def divergence_matrix(self, domain):
-        """The divergence matrix
-
-        Parameters
-        ----------
-        domain : list
-            The domain(s) in which to compute the divergence matrix
-        """
-        raise NotImplementedError
diff --git a/pybamm/discretisations/finite_volume_discretisations.py b/pybamm/discretisations/finite_volume_discretisations.py
index bfa015d6..b97eae23 100644
--- a/pybamm/discretisations/finite_volume_discretisations.py
+++ b/pybamm/discretisations/finite_volume_discretisations.py
@@ -9,33 +9,112 @@ import numpy as np
 from scipy.sparse import spdiags
 
 
-class FiniteVolumeDiscretisation(pybamm.MatrixVectorDiscretisation):
+class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
     """Discretisation using Finite Volumes.
-    Inherits from :class:`pybamm.MatrixVectorDiscretisation`, so we only need to
-    implement the gradient and divergence matrices
 
     Parameters
     ----------
-     mesh : :class:`BaseMesh` (or subclass)
-            The underlying mesh for discretisation
+    mesh : :class:`pybamm.BaseMesh` (or subclass)
+        The underlying mesh for discretisation
 
-    **Extends:** :class:`pybamm.MatrixVectorDiscretisation`
+    **Extends:** :class:`pybamm.BaseDiscretisation`
     """
 
     def __init__(self, mesh):
         super().__init__(mesh)
 
+    def gradient(self, symbol, y_slices, boundary_conditions):
+        """Matrix-vector multiplication to implement the gradient operator.
+        See :meth:`pybamm.BaseDiscretisation.gradient`
+        """
+        # Check that boundary condition keys are hashes (ids)
+        for key in boundary_conditions.keys():
+            assert isinstance(key, int), TypeError(
+                "boundary condition keys should be hashes, not {}".format(type(key))
+            )
+        # Discretise symbol
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        domain = symbol.domain
+        # Add Dirichlet boundary conditions, if defined
+        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)
+            domain = (
+                [domain[0] + "_left ghost cell"]
+                + domain
+                + [domain[-1] + "_right ghost cell"]
+            )
+        gradient_matrix = self.gradient_matrix(domain)
+        return gradient_matrix * discretised_symbol
+
+    def add_ghost_nodes(self, discretised_symbol, lbc, rbc):
+        """
+        Add Dirichlet boundary conditions via ghost nodes.
+
+        For a boundary condition "y = a at the left-hand boundary",
+        we concatenate a ghost node to the start of the vector y with value "2*a - y1"
+        where y1 is the value of the first node.
+        Similarly for the right-hand boundary condition.
+
+        Currently, Dirichlet boundary conditions can only be applied on state
+        variables (e.g. concentration, temperature), and not on expressions.
+        To access the value of the first node (y1), we create a "first_node" object
+        which is a StateVector whose y_slice is the start of the y_slice of
+        discretised_symbol.
+        Similarly, the last node is a StateVector whose y_slice is the end of the
+        y_slice of discretised_symbol
+
+        Parameters
+        ----------
+        discretised_symbol : :class:`pybamm.StateVector` (size n)
+            The discretised variable (a state vector) to which to add ghost nodes
+        lbc : :class:`pybamm.Scalar`
+            Dirichlet bouncary condition on the left-hand side
+        rbc : :class:`pybamm.Scalar`
+            Dirichlet bouncary condition on the right-hand side
+
+        Returns
+        -------
+        :class:`pybamm.Concatenation` (size n+2)
+            Concatenation of the variable (a state vector) and ghost nodes
+
+        """
+        assert isinstance(discretised_symbol, pybamm.StateVector), NotImplementedError(
+            """discretised_symbol must be a StateVector, not {}""".format(
+                type(discretised_symbol)
+            )
+        )
+        # left ghost cell
+        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 self.concatenate(left_ghost_cell, discretised_symbol, right_ghost_cell)
+
     def gradient_matrix(self, domain):
         """
         Gradient matrix for finite volumes in the appropriate domain.
         Equivalent to grad(y) = (y[1:] - y[:-1])/dx
-        See :meth:`pybamm.MatrixVectorDiscretisation.gradient_matrix()`
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) in which to compute the gradient matrix
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) finite volume gradient matrix for the domain
         """
-        assert len(domain) == 1
-        # TODO: implement for when there are several domains
+        # Create appropriate submesh by combining submeshes in domain
+        submesh = self.mesh.combine_submeshes(*domain)
 
-        # implementation for a single domain
-        submesh = self.mesh[domain[0]]
+        # Create matrix using submesh
         n = submesh.npts
         e = 1 / submesh.d_nodes
         data = np.vstack(
@@ -45,17 +124,44 @@ class FiniteVolumeDiscretisation(pybamm.MatrixVectorDiscretisation):
         matrix = spdiags(data, diags, n - 1, n)
         return pybamm.Matrix(matrix)
 
+    def divergence(self, symbol, y_slices, boundary_conditions):
+        """Matrix-vector multiplication to implement the divergence operator.
+        See :meth:`pybamm.BaseDiscretisation.gradient`
+        """
+        # Check that boundary condition keys are hashes (ids)
+        for key in boundary_conditions.keys():
+            assert isinstance(key, int), TypeError(
+                "boundary condition keys should be hashes, not {}".format(type(key))
+            )
+        # Discretise symbol
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        # Add Neumann boundary conditions if defined
+        if symbol.id in boundary_conditions:
+            lbc = boundary_conditions[symbol.id]["left"]
+            rbc = boundary_conditions[symbol.id]["right"]
+            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
+        divergence_matrix = self.divergence_matrix(symbol.domain)
+        return divergence_matrix * discretised_symbol
+
     def divergence_matrix(self, domain):
         """
         Divergence matrix for finite volumes in the appropriate domain.
         Equivalent to div(N) = (N[1:] - N[:-1])/dx
-        See :meth:`pybamm.MatrixVectorDiscretisation.divergence_matrix()`
+
+        Parameters
+        ----------
+        domain : list
+            The domain(s) in which to compute the divergence matrix
+
+        Returns
+        -------
+        :class:`pybamm.Matrix`
+            The (sparse) finite volume divergence matrix for the domain
         """
-        assert len(domain) == 1
-        # TODO: implement for when there are several domains
+        # Create appropriate submesh by combining submeshes in domain
+        submesh = self.mesh.combine_submeshes(*domain)
 
-        # implementation for a single domain
-        submesh = self.mesh[domain[0]]
+        # Create matrix using submesh
         n = submesh.npts + 1
         e = 1 / submesh.d_edges
         data = np.vstack(
diff --git a/pybamm/discretisations/finite_volume_meshes.py b/pybamm/discretisations/finite_volume_meshes.py
index 5a36edf2..08ca3091 100644
--- a/pybamm/discretisations/finite_volume_meshes.py
+++ b/pybamm/discretisations/finite_volume_meshes.py
@@ -51,15 +51,65 @@ class FiniteVolumeMacroMesh(pybamm.BaseMesh):
         self.total_mesh_points = (
             self.neg_mesh_points + (self.sep_mesh_points - 2) + self.pos_mesh_points
         )
-        self["whole cell"] = self.submeshclass(
-            np.concatenate(
-                [
-                    self["negative electrode"].edges,
-                    self["separator"].edges[1:-1],
-                    self["positive electrode"].edges,
-                ]
-            )
+        self["whole cell"] = self.combine_submeshes(
+            "negative electrode", "separator", "positive electrode"
+        )
+
+        # Add ghost meshes for ghost cells for Dirichlet boundary conditions
+        self.add_ghost_meshes()
+
+    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
+        not aligned).
+
+        Parameters
+        ----------
+        submeshnames: list of str
+            The names of the submeshes to be combined
+
+        Returns
+        -------
+        submesh: :class:`self.submeshclass`
+            A new submesh with the class defined by self.submeshclass
+        """
+        # 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
+            combined_submesh_edges = np.concatenate(
+                [self[submeshnames[0]].edges]
+                + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]]
+            )
+            return self.submeshclass(combined_submesh_edges)
+        else:
+            raise pybamm.DomainError("submesh edges are not aligned")
+
+    def add_ghost_meshes(self):
+        """
+        Create meshes for potential ghost nodes on either side of each submesh, using
+        self.submeshclass
+        This will be useful for calculating the gradient with Dirichlet BCs.
+        """
+        # Get all submeshes relating to space (i.e. exclude time)
+        submeshes = [
+            (name, submesh) for name, submesh in self.items() if name != "time"
+        ]
+        for submeshname, 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[submeshname + "_left ghost cell"] = self.submeshclass(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"] = self.submeshclass(rgs_edges)
 
 
 class FiniteVolumeSubmesh:

discritisation of boundary conditions
I was just looking at the discretisation of the boundary conditions. At the moment this seems to be mainly done in the `gradient` and `divergence` functions of `MatrixVectorDiscretisation`, which both do this:

```python
       # Add boundary conditions if defined
        if symbol.id in boundary_conditions:
            lbc, rbc = boundary_conditions[symbol.id]
            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
        gradient_matrix = self.gradient_matrix(domain)
        return gradient_matrix * discretised_symbol
```
So it looks like extra ghost nodes are added to the current state vector to implement the bc. Should this be done instead when the `pybamm.Variable` node is processed in `process_symbol`? Also, the way the boundary conditions are implemented would depend on the method, so this should probably be done in `FiniteVolumeDiscretisation`? @tinosulzer? 


diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index adfe5b30..7ff51bea 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -15,19 +15,23 @@ class MeshForTesting(pybamm.BaseMesh):
         self["positive electrode"] = self.submeshclass(self["whole cell"].nodes[40:])
 
 
-class DiscretisationForTesting(pybamm.MatrixVectorDiscretisation):
-    """Interpolating operators."""
+class DiscretisationForTesting(pybamm.BaseDiscretisation):
+    """Identity operators, no boundary conditions."""
 
     def __init__(self, mesh):
         super().__init__(mesh)
 
-    def gradient_matrix(self, domain):
-        n = self.mesh[domain[0]].npts
-        return pybamm.Matrix(np.eye(n))
+    def gradient(self, symbol, y_slices, boundary_conditions):
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        n = self.mesh[symbol.domain[0]].npts
+        gradient_matrix = pybamm.Matrix(np.eye(n))
+        return gradient_matrix * discretised_symbol
 
-    def divergence_matrix(self, domain):
-        n = self.mesh[domain[0]].npts
-        return pybamm.Matrix(np.eye(n))
+    def divergence(self, symbol, y_slices, boundary_conditions):
+        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+        n = self.mesh[symbol.domain[0]].npts
+        divergence_matrix = pybamm.Matrix(np.eye(n))
+        return divergence_matrix * discretised_symbol
 
 
 class TestDiscretise(unittest.TestCase):
@@ -165,18 +169,13 @@ class TestDiscretise(unittest.TestCase):
             disc.gradient(None, None, {})
         with self.assertRaises(NotImplementedError):
             disc.divergence(None, None, {})
-        disc = pybamm.MatrixVectorDiscretisation(None)
-        with self.assertRaises(NotImplementedError):
-            disc.gradient_matrix(None)
-        with self.assertRaises(NotImplementedError):
-            disc.divergence_matrix(None)
 
     def test_process_initial_conditions(self):
         # one equation
         c = pybamm.Variable("c", domain=["whole cell"])
         initial_conditions = {c: pybamm.Scalar(3)}
         mesh = MeshForTesting()
-        disc = DiscretisationForTesting(mesh)
+        disc = pybamm.BaseDiscretisation(mesh)
         y0 = disc.process_initial_conditions(initial_conditions)
         np.testing.assert_array_equal(
             y0[c].evaluate(0, None), 3 * np.ones_like(mesh["whole cell"].nodes)
diff --git a/tests/test_discretisations/test_finite_volume_discretisations.py b/tests/test_discretisations/test_finite_volume_discretisations.py
index 1d00b816..f0c444ed 100644
--- a/tests/test_discretisations/test_finite_volume_discretisations.py
+++ b/tests/test_discretisations/test_finite_volume_discretisations.py
@@ -8,7 +8,92 @@ import unittest
 
 
 class TestFiniteVolumeDiscretisation(unittest.TestCase):
-    def test_grad_div_shapes(self):
+    def test_add_ghost_nodes(self):
+        # Set up
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+        # Add ghost nodes
+        var = pybamm.Variable("var", domain=["whole cell"])
+        y_slices = disc.get_variable_slices([var])
+        discretised_symbol = pybamm.StateVector(y_slices[var.id])
+        lbc = pybamm.Scalar(0)
+        rbc = pybamm.Scalar(3)
+        symbol_plus_ghost = disc.add_ghost_nodes(discretised_symbol, lbc, rbc)
+
+        # Test
+        y_test = np.ones_like(mesh["whole cell"].nodes)
+        np.testing.assert_array_equal(
+            symbol_plus_ghost.evaluate(None, y_test)[1:-1],
+            discretised_symbol.evaluate(None, y_test),
+        )
+        self.assertEqual(
+            (
+                symbol_plus_ghost.evaluate(None, y_test)[0]
+                + symbol_plus_ghost.evaluate(None, y_test)[1]
+            )
+            / 2,
+            0,
+        )
+        self.assertEqual(
+            (
+                symbol_plus_ghost.evaluate(None, y_test)[-2]
+                + symbol_plus_ghost.evaluate(None, y_test)[-1]
+            )
+            / 2,
+            3,
+        )
+
+    def test_grad_div_shapes_Dirichlet_bcs(self):
+        """
+        Test grad and div with Dirichlet boundary conditions (applied by grad on var)
+        """
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+        # grad
+        var = pybamm.Variable("var", domain=["whole cell"])
+        grad_eqn = pybamm.grad(var)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
+        }
+        y_slices = disc.get_variable_slices([var])
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+
+        constant_y = np.ones_like(mesh["whole cell"].nodes)
+        np.testing.assert_array_equal(
+            grad_eqn_disc.evaluate(None, constant_y),
+            np.zeros_like(mesh["whole cell"].edges),
+        )
+
+        # div: test on linear y (should have laplacian zero) so change bcs
+        linear_y = mesh["whole cell"].nodes
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}
+        }
+
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            grad_eqn_disc.evaluate(None, linear_y),
+            np.ones_like(mesh["whole cell"].edges),
+        )
+
+        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(None, linear_y),
+            np.zeros_like(mesh["whole cell"].nodes),
+        )
+
+    def test_grad_div_shapes_Neumann_bcs(self):
+        """Test grad and div with Neumann boundary conditions (applied by div on N)"""
         param = pybamm.ParameterValues(
             base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
         )
@@ -46,7 +131,54 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             np.zeros_like(mesh["whole cell"].nodes),
         )
 
-    def test_grad_convergence(self):
+    def test_grad_div_shapes_mixed_domain(self):
+        """
+        Test grad and div with Dirichlet boundary conditions (applied by grad on var)
+        """
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 2)
+        disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+        # grad
+        var = pybamm.Variable("var", domain=["negative electrode", "separator"])
+        grad_eqn = pybamm.grad(var)
+        boundary_conditions = {
+            var.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
+        }
+        y_slices = disc.get_variable_slices([var])
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+
+        combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
+        constant_y = np.ones_like(combined_submesh.nodes)
+        np.testing.assert_array_equal(
+            grad_eqn_disc.evaluate(None, constant_y),
+            np.zeros_like(combined_submesh.edges),
+        )
+
+        # div: test on linear y (should have laplacian zero) so change bcs
+        linear_y = combined_submesh.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]),
+            }
+        }
+
+        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            grad_eqn_disc.evaluate(None, linear_y), np.ones_like(combined_submesh.edges)
+        )
+
+        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        np.testing.assert_array_almost_equal(
+            div_eqn_disc.evaluate(None, linear_y), np.zeros_like(combined_submesh.nodes)
+        )
+
+    def test_grad_convergence_without_bcs(self):
         # Convergence
         param = pybamm.ParameterValues(
             base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
@@ -54,10 +186,8 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         var = pybamm.Variable("var", domain=["whole cell"])
         grad_eqn = pybamm.grad(var)
 
-        # Prepare convergence testing
-        ns = [50, 100, 200]
-        errs = [0] * len(ns)
-        for i, n in enumerate(ns):
+        # Function for convergence testing
+        def get_l2_error(n):
             # Set up discretisation
             mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
             disc = pybamm.FiniteVolumeDiscretisation(mesh)
@@ -72,14 +202,96 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             grad_approx = grad_eqn_disc.evaluate(None, y)
 
             # Calculate errors
-            errs[i] = np.linalg.norm(grad_approx - grad_exact) / np.linalg.norm(
-                grad_exact
-            )
+            return np.linalg.norm(grad_approx - grad_exact) / np.linalg.norm(grad_exact)
 
-        # Expect h**2 convergence
-        [self.assertLess(errs[i + 1] / errs[i], 0.26) for i in range(len(errs) - 1)]
+        # Get errors
+        ns = 100 * (2 ** np.arange(1, 7))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
+
+        # Get rates: expect h**2 convergence
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
+
+    def test_grad_convergence_with_bcs(self):
+        # Convergence
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        var = pybamm.Variable("var", domain=["whole cell"])
+        grad_eqn = pybamm.grad(var)
+        boundary_conditions = {
+            var.id: {
+                "left": pybamm.Scalar(np.sin(0)),
+                "right": pybamm.Scalar(np.sin(1)),
+            }
+        }
+
+        # Function for convergence testing
+        def get_l2_error(n):
+            # Set up discretisation
+            mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
+            disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+            # Define exact solutions
+            y = np.sin(mesh["whole cell"].nodes)
+            grad_exact = np.cos(mesh["whole cell"].edges)
+
+            # Discretise and evaluate
+            y_slices = disc.get_variable_slices([var])
+            grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+            grad_approx = grad_eqn_disc.evaluate(None, y)
+
+            # Calculate errors
+            return np.linalg.norm(grad_approx - grad_exact) / np.linalg.norm(grad_exact)
+
+        # Get errors
+        ns = 100 * (2 ** np.arange(1, 7))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
+
+        # Get rates: expect h**1.5 convergence
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.49 * np.ones_like(rates), rates)
+
+    def test_div_convergence_internal(self):
+        # Convergence
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.1, "Ls": 0.2, "Lp": 0.3}
+        )
+        var = pybamm.Variable("var", domain=["whole cell"])
+        N = pybamm.grad(var)
+        div_eqn = pybamm.div(N)
+        boundary_conditions = {
+            N.id: {"left": pybamm.Scalar(np.cos(0)), "right": pybamm.Scalar(np.cos(1))}
+        }
+
+        # Function for convergence testing
+        def get_l2_error(n):
+            # Set up discretisation
+            mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
+            disc = pybamm.FiniteVolumeDiscretisation(mesh)
+
+            # Define exact solutions
+            y = np.sin(mesh["whole cell"].nodes)
+            div_exact_internal = -np.sin(mesh["whole cell"].nodes[1:-1])
+
+            # Discretise and evaluate
+            y_slices = disc.get_variable_slices([var])
+            div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+            div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
+
+            # Calculate errors
+            return np.linalg.norm(
+                div_approx_internal - div_exact_internal
+            ) / np.linalg.norm(div_exact_internal)
+
+        # Get errors
+        ns = 10 * (2 ** np.arange(1, 6))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
+
+        # Get rates: expect h**2 convergence
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
 
-    @unittest.skip("div errors do not converge as expected")
     def test_div_convergence(self):
         # Convergence
         param = pybamm.ParameterValues(
@@ -92,28 +304,33 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             N.id: {"left": pybamm.Scalar(np.cos(0)), "right": pybamm.Scalar(np.cos(1))}
         }
 
-        # Prepare convergence testing
-        ns = [50, 100, 200]
-        errs = [0] * len(ns)
-        for i, n in enumerate(ns):
+        # Function for convergence testing
+        def get_l2_error(n):
             # Set up discretisation
             mesh = pybamm.FiniteVolumeMacroMesh(param, target_npts=n)
             disc = pybamm.FiniteVolumeDiscretisation(mesh)
 
             # Define exact solutions
             y = np.sin(mesh["whole cell"].nodes)
-            div_exact = -np.sin(mesh["whole cell"].nodes)
+            div_exact_internal = -np.sin(mesh["whole cell"].nodes[1:-1])
 
             # Discretise and evaluate
             y_slices = disc.get_variable_slices([var])
             div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
-            div_approx = div_eqn_disc.evaluate(None, y)
+            div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
 
             # Calculate errors
-            errs[i] = np.linalg.norm(div_approx - div_exact) / np.linalg.norm(div_exact)
+            return np.linalg.norm(
+                div_approx_internal - div_exact_internal
+            ) / np.linalg.norm(div_exact_internal)
+
+        # Get errors
+        ns = 10 * (2 ** np.arange(1, 6))
+        errs = np.array([get_l2_error(int(n)) for n in ns])
 
-        # Expect h**2 convergence
-        [self.assertLess(errs[i + 1] / errs[i], 0.26) for i in range(len(errs) - 1)]
+        # Get rates: expect h**1.5 convergence because of boundary conditions
+        rates = np.log2(errs[:-1] / errs[1:])
+        np.testing.assert_array_less(1.49 * np.ones_like(rates), rates)
 
 
 if __name__ == "__main__":
diff --git a/tests/test_discretisations/test_finite_volume_meshes.py b/tests/test_discretisations/test_finite_volume_meshes.py
index 983a83b4..42284f68 100644
--- a/tests/test_discretisations/test_finite_volume_meshes.py
+++ b/tests/test_discretisations/test_finite_volume_meshes.py
@@ -50,6 +50,52 @@ class TestFiniteVolumeMesh(unittest.TestCase):
         self.assertEqual(len(mesh["positive electrode"].edges), mesh.pos_mesh_points)
         self.assertEqual(len(mesh["whole cell"].edges), mesh.total_mesh_points)
 
+    def test_combine_submeshes(self):
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.01, "Ls": 0.5, "Lp": 0.12}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 50)
+        submesh = mesh.combine_submeshes("negative electrode", "separator")
+        self.assertEqual(submesh.edges[0], 0)
+        self.assertEqual(submesh.edges[-1], mesh["separator"].edges[-1])
+        self.assertAlmostEqual(
+            np.linalg.norm(
+                submesh.nodes
+                - np.concatenate(
+                    [mesh["negative electrode"].nodes, mesh["separator"].nodes]
+                )
+            ),
+            0,
+        )
+        with self.assertRaises(pybamm.DomainError):
+            submesh = mesh.combine_submeshes("negative electrode", "positive electrode")
+
+    def test_ghost_cells(self):
+        param = pybamm.ParameterValues(
+            base_parameters={"Ln": 0.01, "Ls": 0.5, "Lp": 0.12}
+        )
+        mesh = pybamm.FiniteVolumeMacroMesh(param, 50)
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges,
+            mesh["whole cell_left ghost cell"].edges,
+        )
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges[1],
+            mesh["negative electrode"].edges[0],
+        )
+        np.testing.assert_array_equal(
+            mesh["negative electrode_left ghost cell"].edges[0],
+            -mesh["negative electrode"].edges[1],
+        )
+        np.testing.assert_array_equal(
+            mesh["positive electrode_right ghost cell"].edges,
+            mesh["whole cell_right ghost cell"].edges,
+        )
+        np.testing.assert_array_equal(
+            mesh["positive electrode_right ghost cell"].edges[0],
+            mesh["positive electrode"].edges[-1],
+        )
+
 
 if __name__ == "__main__":
     print("Add -v for more debug output")

pip install -e .[dev,docs]
parse_log_pytest
requirements.txt
pytest
3.6
requirements.txt requirements-dev.txt
pytest --no-header -rA --tb=line --color=no -p no:cacheprovider -W ignore::DeprecationWarning
alabaster==0.7.13
anyio==3.6.2
anytree==2.8.0
argon2-cffi==21.3.0
argon2-cffi-bindings==21.2.0
async-generator==1.10
attrs==22.2.0
Babel==2.11.0
backcall==0.2.0
bleach==4.1.0
certifi==2021.5.30
cffi==1.15.1
charset-normalizer==2.0.12
comm==0.1.4
contextvars==2.4
cycler==0.11.0
dataclasses==0.8
decorator==5.1.1
defusedxml==0.7.1
docutils==0.18.1
entrypoints==0.4
flake8==3.9.2
guzzle-sphinx-theme==0.7.11
idna==3.10
imagesize==1.4.1
immutables==0.19
importlib-metadata==4.8.3
iniconfig==1.1.1
ipykernel==5.5.6
ipython==7.16.3
ipython-genutils==0.2.0
ipywidgets==7.8.5
jedi==0.17.2
Jinja2==3.0.3
json5==0.9.16
jsonschema==3.2.0
jupyter==1.1.1
jupyter-client==7.1.2
jupyter-console==6.4.3
jupyter-core==4.9.2
jupyter-server==1.13.1
jupyterlab==3.2.9
jupyterlab-pygments==0.1.2
jupyterlab-server==2.10.3
jupyterlab_widgets==1.1.11
kiwisolver==1.3.1
MarkupSafe==2.0.1
matplotlib==3.3.4
mccabe==0.6.1
mistune==0.8.4
nbclassic==0.3.5
nbclient==0.5.9
nbconvert==6.0.7
nbformat==5.1.3
nest-asyncio==1.6.0
notebook==6.4.10
numpy==1.19.5
packaging==21.3
pandas==1.1.5
pandocfilters==1.5.1
parso==0.7.1
pexpect==4.9.0
pickleshare==0.7.5
Pillow==8.4.0
pluggy==1.0.0
prometheus-client==0.17.1
prompt-toolkit==3.0.36
ptyprocess==0.7.0
py==1.11.0
-e git+https://github.com/tinosulzer/PyBaMM.git@a10e4d5a937cf685c02c7c2f40c52c43b3fd51e8#egg=pybamm
pycodestyle==2.7.0
pycparser==2.21
pyflakes==2.3.1
Pygments==2.14.0
pyparsing==3.1.4
pyrsistent==0.18.0
pytest==7.0.1
python-dateutil==2.9.0.post0
pytz==2025.2
pyzmq==25.1.2
requests==2.27.1
scipy==1.5.4
Send2Trash==1.8.3
six==1.17.0
sniffio==1.2.0
snowballstemmer==2.2.0
Sphinx==5.3.0
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
terminado==0.12.1
testpath==0.6.0
tomli==1.2.3
tornado==6.1
traitlets==4.3.3
typing_extensions==4.1.1
urllib3==1.26.20
wcwidth==0.2.13
webencodings==0.5.1
websocket-client==1.3.1
widgetsnbextension==3.6.10
zipp==3.6.0

name: PyBaMM
channels:
  - defaults
  - https://repo.anaconda.com/pkgs/main
  - https://repo.anaconda.com/pkgs/r
  - conda-forge
dependencies:
  - _libgcc_mutex=0.1=main
  - _openmp_mutex=5.1=1_gnu
  - ca-certificates=2025.2.25=h06a4308_0
  - certifi=2021.5.30=py36h06a4308_0
  - ld_impl_linux-64=2.40=h12ee557_0
  - libffi=3.3=he6710b0_2
  - libgcc-ng=11.2.0=h1234567_1
  - libgomp=11.2.0=h1234567_1
  - libstdcxx-ng=11.2.0=h1234567_1
  - ncurses=6.4=h6a678d5_0
  - openssl=1.1.1w=h7f8727e_0
  - pip=21.2.2=py36h06a4308_0
  - python=3.6.13=h12debd9_1
  - readline=8.2=h5eee18b_0
  - setuptools=58.0.4=py36h06a4308_0
  - sqlite=3.45.3=h5eee18b_0
  - tk=8.6.14=h39e8969_0
  - wheel=0.37.1=pyhd3eb1b0_0
  - xz=5.6.4=h5eee18b_1
  - zlib=1.2.13=h5eee18b_1
  - pip:
      - alabaster==0.7.13
      - anyio==3.6.2
      - anytree==2.8.0
      - argon2-cffi==21.3.0
      - argon2-cffi-bindings==21.2.0
      - async-generator==1.10
      - attrs==22.2.0
      - babel==2.11.0
      - backcall==0.2.0
      - bleach==4.1.0
      - cffi==1.15.1
      - charset-normalizer==2.0.12
      - comm==0.1.4
      - contextvars==2.4
      - cycler==0.11.0
      - dataclasses==0.8
      - decorator==5.1.1
      - defusedxml==0.7.1
      - docutils==0.18.1
      - entrypoints==0.4
      - flake8==3.9.2
      - guzzle-sphinx-theme==0.7.11
      - idna==3.10
      - imagesize==1.4.1
      - immutables==0.19
      - importlib-metadata==4.8.3
      - iniconfig==1.1.1
      - ipykernel==5.5.6
      - ipython==7.16.3
      - ipython-genutils==0.2.0
      - ipywidgets==7.8.5
      - jedi==0.17.2
      - jinja2==3.0.3
      - json5==0.9.16
      - jsonschema==3.2.0
      - jupyter==1.1.1
      - jupyter-client==7.1.2
      - jupyter-console==6.4.3
      - jupyter-core==4.9.2
      - jupyter-server==1.13.1
      - jupyterlab==3.2.9
      - jupyterlab-pygments==0.1.2
      - jupyterlab-server==2.10.3
      - jupyterlab-widgets==1.1.11
      - kiwisolver==1.3.1
      - markupsafe==2.0.1
      - matplotlib==3.3.4
      - mccabe==0.6.1
      - mistune==0.8.4
      - nbclassic==0.3.5
      - nbclient==0.5.9
      - nbconvert==6.0.7
      - nbformat==5.1.3
      - nest-asyncio==1.6.0
      - notebook==6.4.10
      - numpy==1.19.5
      - packaging==21.3
      - pandas==1.1.5
      - pandocfilters==1.5.1
      - parso==0.7.1
      - pexpect==4.9.0
      - pickleshare==0.7.5
      - pillow==8.4.0
      - pluggy==1.0.0
      - prometheus-client==0.17.1
      - prompt-toolkit==3.0.36
      - ptyprocess==0.7.0
      - py==1.11.0
      - pycodestyle==2.7.0
      - pycparser==2.21
      - pyflakes==2.3.1
      - pygments==2.14.0
      - pyparsing==3.1.4
      - pyrsistent==0.18.0
      - pytest==7.0.1
      - python-dateutil==2.9.0.post0
      - pytz==2025.2
      - pyzmq==25.1.2
      - requests==2.27.1
      - scipy==1.5.4
      - send2trash==1.8.3
      - six==1.17.0
      - sniffio==1.2.0
      - snowballstemmer==2.2.0
      - sphinx==5.3.0
      - sphinxcontrib-applehelp==1.0.2
      - sphinxcontrib-devhelp==1.0.2
      - sphinxcontrib-htmlhelp==2.0.0
      - sphinxcontrib-jsmath==1.0.1
      - sphinxcontrib-qthelp==1.0.3
      - sphinxcontrib-serializinghtml==1.1.5
      - terminado==0.12.1
      - testpath==0.6.0
      - tomli==1.2.3
      - tornado==6.1
      - traitlets==4.3.3
      - typing-extensions==4.1.1
      - urllib3==1.26.20
      - wcwidth==0.2.13
      - webencodings==0.5.1
      - websocket-client==1.3.1
      - widgetsnbextension==3.6.10
      - zipp==3.6.0
prefix: /opt/conda/envs/PyBaMM

tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_add_ghost_nodes tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_grad_convergence_with_bcs tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_grad_div_shapes_Dirichlet_bcs tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_grad_div_shapes_mixed_domain tests/test_discretisations/test_finite_volume_meshes.py::TestFiniteVolumeMesh::test_combine_submeshes tests/test_discretisations/test_finite_volume_meshes.py::TestFiniteVolumeMesh::test_ghost_cells
tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_concatenation tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_concatenation_of_scalars tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_core_NotImplementedErrors tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_discretise_slicing tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_discretise_spatial_operator tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_complex_expression tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_dict tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_initial_conditions tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_model tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_process_symbol_base tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_scalar_to_vector tests/test_discretisations/test_base_discretisations.py::TestDiscretise::test_vector_of_ones tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_div_convergence tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_div_convergence_internal tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_grad_convergence_without_bcs tests/test_discretisations/test_finite_volume_discretisations.py::TestFiniteVolumeDiscretisation::test_grad_div_shapes_Neumann_bcs tests/test_discretisations/test_finite_volume_meshes.py::TestFiniteVolumeMesh::test_mesh_creation tests/test_discretisations/test_finite_volume_meshes.py::TestFiniteVolumeMesh::test_mesh_sizes
docs/source/discretisations/core.rst pybamm/__init__.py docs/source/discretisations/base_discretisation.rst .flake8 pybamm/discretisations/finite_volume_discretisations.py .gitignore pybamm/discretisations/finite_volume_meshes.py docs/source/discretisations/index.rst pybamm/discretisations/base_discretisation.py
docs/source/discretisations/core.rst pybamm/__init__.py docs/source/discretisations/base_discretisation.rst .flake8 pybamm/discretisations/finite_volume_discretisations.py .gitignore pybamm/discretisations/finite_volume_meshes.py docs/source/discretisations/index.rst pybamm/discretisations/base_discretisation.py