Repository: tinosulzer/PyBaMM

Issue:
enable per-equation discretiation
depends on #128 

dict of {variable: discretisation} passed as input to base discretisation class. 

The discretiation class passed in the dict handles things like gradient/divergence discretisation + broadcasting

Discussion:

Patch:
diff --git a/docs/index.rst b/docs/index.rst
index 216607ff..a145fadb 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -28,4 +28,6 @@ Contents
     source/models/index
     source/parameters/index
     source/discretisations/index
+    source/spatial_methods/index
     source/solvers/index
+
diff --git a/docs/source/discretisations/base_discretisation.rst b/docs/source/discretisations/base_discretisation.rst
deleted file mode 100644
index 95531920..00000000
--- a/docs/source/discretisations/base_discretisation.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-Base Discretisation
-===================
-
-.. autoclass:: pybamm.BaseDiscretisation
-  :members:
diff --git a/docs/source/discretisations/discretisation.rst b/docs/source/discretisations/discretisation.rst
new file mode 100644
index 00000000..af2e06f2
--- /dev/null
+++ b/docs/source/discretisations/discretisation.rst
@@ -0,0 +1,5 @@
+Discretisation
+==============
+
+.. autoclass:: pybamm.Discretisation
+  :members:
diff --git a/docs/source/discretisations/finite_volume_discretisations.rst b/docs/source/discretisations/finite_volume_discretisations.rst
deleted file mode 100644
index 3d16c187..00000000
--- a/docs/source/discretisations/finite_volume_discretisations.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-Finite Volume Discretisations
-=============================
-
-.. autoclass:: pybamm.FiniteVolumeDiscretisation
-  :members:
-
-.. autoclass:: pybamm.NodeToEdge
-  :members:
diff --git a/docs/source/discretisations/index.rst b/docs/source/discretisations/index.rst
index 4ad52af9..3f9c1846 100644
--- a/docs/source/discretisations/index.rst
+++ b/docs/source/discretisations/index.rst
@@ -3,5 +3,4 @@ Discretisations
 
 .. toctree::
 
-  base_discretisation
-  finite_volume_discretisations
+  discretisation
diff --git a/docs/source/spatial_methods/finite_volume.rst b/docs/source/spatial_methods/finite_volume.rst
new file mode 100644
index 00000000..f57c5e9e
--- /dev/null
+++ b/docs/source/spatial_methods/finite_volume.rst
@@ -0,0 +1,6 @@
+Finite Volume
+=============
+
+.. autoclass:: pybamm.FiniteVolume
+    :members:
+
diff --git a/docs/source/spatial_methods/index.rst b/docs/source/spatial_methods/index.rst
new file mode 100644
index 00000000..f83702cc
--- /dev/null
+++ b/docs/source/spatial_methods/index.rst
@@ -0,0 +1,7 @@
+Spatial Methods
+===============
+
+.. toctree::
+
+  spatial_method
+  finite_volume
diff --git a/docs/source/spatial_methods/spatial_method.rst b/docs/source/spatial_methods/spatial_method.rst
new file mode 100644
index 00000000..72d243dd
--- /dev/null
+++ b/docs/source/spatial_methods/spatial_method.rst
@@ -0,0 +1,5 @@
+Spatial Method
+==============
+
+.. autoclass:: pybamm.SpatialMethod
+    :members:
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 7d9c9a77..ac4b6da7 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -124,14 +124,17 @@ from .geometry.geometry import (
 #
 # Mesh and Discretisation classes
 #
-from .discretisations.base_discretisation import BaseDiscretisation
-from .discretisations.finite_volume_discretisations import (
-    FiniteVolumeDiscretisation,
-    NodeToEdge,
-)
+from .discretisations.discretisation import Discretisation
+from .meshes.meshes import KNOWN_DOMAINS
 from .meshes.meshes import Mesh
 from .meshes.submeshes import SubMesh1D, Uniform1DSubMesh
 
+#
+# Spatial Methods
+#
+from .spatial_methods.spatial_method import SpatialMethod
+from .spatial_methods.finite_volume import FiniteVolume, NodeToEdge
+
 #
 # Simulation class
 #
diff --git a/pybamm/discretisations/base_discretisation.py b/pybamm/discretisations/discretisation.py
similarity index 61%
rename from pybamm/discretisations/base_discretisation.py
rename to pybamm/discretisations/discretisation.py
index 067cfdb7..a572b4fb 100644
--- a/pybamm/discretisations/base_discretisation.py
+++ b/pybamm/discretisations/discretisation.py
@@ -9,22 +9,27 @@ import copy
 import numpy as np
 
 
-class BaseDiscretisation(object):
-    """The base discretisation class, with methods to process a model and replace
+class Discretisation(object):
+    """The discretisation class, with methods to process a model and replace
     Spatial Operators with Matrices and Variables with StateVectors
 
     Parameters
     ----------
-    mesh_type: class name
-        the type of combined mesh being used (e.g. Mesh)
-    submesh_pts : dict
-        the number of points on each of the subdomains
-    submesh_types : dict
-        the types of submeshes being employed on each of the subdomains
+    mesh : pybamm.Mesh
+            contains all submeshes to be used on each domain
+    spatial_methods : dict
+            a dictionary of the spatial method to be used on each
+            domain. The keys correspond to the keys in a pybamm.Model
     """
 
-    def __init__(self, mesh):
+    def __init__(self, mesh, spatial_methods):
         self._mesh = mesh
+        self._spatial_methods = {
+            dom: method(mesh) for dom, method in spatial_methods.items()
+        }
+        self._bcs = {}
+        self._y_slices = {}
+        self._variables = []
 
     @property
     def mesh(self):
@@ -41,44 +46,35 @@ class BaseDiscretisation(object):
             boundary_conditions (all dicts of {variable: equation})
 
         """
+        # set boundary conditions (only need key ids for boundary_conditions)
+        self._bcs = {key.id: value for key, value in model.boundary_conditions.items()}
+        # set variables (we require the full variable not just id)
+        self._variables = list(model.rhs.keys()) + list(model.algebraic.keys())
 
         # Set the y split for variables
-        variables = list(model.rhs.keys()) + list(model.algebraic.keys())
-        y_slices = self.get_variable_slices(variables)
+        self.set_variable_slices()
 
         # Process initial condtions
-        self.process_initial_conditions(model, y_slices)
+        self.process_initial_conditions(model)
 
-        # Process differential and algebraic equations
-        self.process_rhs_and_algebraic(model, y_slices)
+        # Process parabolic and elliptic equations
+        self.process_rhs_and_algebraic(model)
 
         # Discretise variables (applying boundary conditions)
         # Note that we **do not** discretise the keys of model.rhs,
         # model.initial_conditions and model.boundary_conditions
-        model.variables = self.process_dict(
-            model.variables, y_slices, model.boundary_conditions
-        )
+        model.variables = self.process_dict(model.variables)
 
         # Check that resulting model makes sense
         self.check_model(model)
 
-    def get_variable_slices(self, variables):
-        """Set the slicing for variables.
-
-        Parameters
-        ----------
-        variables : list of Variables
-            Variables for which to set the slices
-
-        Returns
-        -------
-        y_slices : dict of {variable id: slice}
-            The slices to take when solving (assigning chunks of y to each vector)
+    def set_variable_slices(self):
+        """Sets the slicing for variables.
         """
-        y_slices = {variable.id: None for variable in variables}
+        y_slices = {variable.id: None for variable in self._variables}
         start = 0
         end = 0
-        for variable in variables:
+        for variable in self._variables:
             # If domain is empty then variable has size 1
             if variable.domain == []:
                 end += 1
@@ -88,10 +84,13 @@ class BaseDiscretisation(object):
                     end += self.mesh[dom].npts
             y_slices[variable.id] = slice(start, end)
             start = end
+        self._y_slices = y_slices
 
-        return y_slices
+        assert isinstance(self._y_slices, dict), ValueError(
+            """y_slices should be dict, not {}""".format(type(self._y_slices))
+        )
 
-    def process_initial_conditions(self, model, y_slices):
+    def process_initial_conditions(self, model):
         """Discretise model initial_conditions.
         Currently inplace, could be changed to return a new model.
 
@@ -100,8 +99,6 @@ class BaseDiscretisation(object):
         model : :class:`pybamm.BaseModel` (or subclass)
             Model to dicretise. Must have attributes rhs, initial_conditions and
             boundary_conditions (all dicts of {variable: equation})
-        y_slices : dict of {variable id: slice}
-            The slices to assign to StateVectors when discretising
 
         """
         # Discretise initial conditions
@@ -109,21 +106,20 @@ class BaseDiscretisation(object):
         model.initial_conditions_ydot = self.process_dict(model.initial_conditions_ydot)
 
         # Concatenate initial conditions into a single vector
-
         # check that all initial conditions are set
         model.concatenated_initial_conditions = self._concatenate_init(
-            model.initial_conditions, y_slices
+            model.initial_conditions
         ).evaluate(0, None)
 
         # evaluate initial conditions for ydot if they exist
         if len(model.initial_conditions_ydot) > 0:
             model.concatenated_initial_conditions_ydot = self._concatenate_init(
-                model.initial_conditions_ydot, y_slices
+                model.initial_conditions_ydot
             ).evaluate(0, None)
         else:
             model.concatenated_initial_conditions_ydot = np.array([])
 
-    def process_rhs_and_algebraic(self, model, y_slices):
+    def process_rhs_and_algebraic(self, model):
         """Discretise model equations - differential ('rhs') and algebraic.
         Currently inplace, could be changed to return a new model.
 
@@ -132,22 +128,17 @@ class BaseDiscretisation(object):
         model : :class:`pybamm.BaseModel` (or subclass)
             Model to dicretise. Must have attributes rhs, initial_conditions and
             boundary_conditions (all dicts of {variable: equation})
-        y_slices : dict of {variable id: slice}
-            The slices to assign to StateVectors when discretising
-
         """
         # Discretise right-hand sides, passing domain from variable
-        model.rhs = self.process_dict(model.rhs, y_slices, model.boundary_conditions)
+        model.rhs = self.process_dict(model.rhs)
         # Concatenate rhs into a single state vector
         model.concatenated_rhs = self.concatenate(*model.rhs.values())
 
         # Discretise and concatenate algebraic equations
-        model.algebraic = self.process_dict(
-            model.algebraic, y_slices, model.boundary_conditions
-        )
+        model.algebraic = self.process_dict(model.algebraic)
         model.concatenated_algebraic = self.concatenate(*model.algebraic.values())
 
-    def process_dict(self, var_eqn_dict, y_slices=None, boundary_conditions={}):
+    def process_dict(self, var_eqn_dict):
         """Discretise a dictionary of {variable: equation}, broadcasting if necessary
         (can be model.rhs, model.initial_conditions or model.variables).
 
@@ -156,11 +147,6 @@ class BaseDiscretisation(object):
         var_eqn_dict : dict
             Equations ({variable: equation} dict) to dicretise
             (can be model.rhs, model.initial_conditions or model.variables)
-        y_slices : dict of {variable id: slice}
-            The slices to assign to StateVectors when discretising
-        boundary_conditions : dict
-            Boundary conditions ({symbol.id: {"left": left bc, "right": right bc}} dict)
-            associated with var_eqn_dict to dicretise
 
         Returns
         -------
@@ -168,34 +154,30 @@ class BaseDiscretisation(object):
             Discretised equations
 
         """
-        # reformat boundary conditions (for rhs and variables)
-        boundary_conditions = {
-            key.id: value for key, value in boundary_conditions.items()
-        }
-        for variable, equation in var_eqn_dict.items():
-            if equation.evaluates_to_number():
-                # Broadcast scalar equation to the domain specified by variable.domain
-                equation = self.broadcast(equation, variable.domain)
 
-            # Process symbol (original or broadcasted)
-            var_eqn_dict[variable] = self.process_symbol(
-                equation, y_slices, boundary_conditions
-            )
+        for eqn_key, eqn in var_eqn_dict.items():
+
+            if eqn.evaluates_to_number():
+                if eqn_key.domain == []:
+                    eqn = pybamm.NumpyBroadcast(eqn, eqn_key.domain, {})
+                else:
+                    eqn = self._spatial_methods[eqn_key.domain[0]].broadcast(
+                        eqn, eqn_key.domain
+                    )
 
+            # Process symbol (original or broadcasted)
+            var_eqn_dict[eqn_key] = self.process_symbol(eqn)
+            # note we are sending in the key.id here so we don't have to
+            # keep calling .id
         return var_eqn_dict
 
-    def process_symbol(self, symbol, y_slices=None, boundary_conditions={}):
+    def process_symbol(self, symbol):
         """Discretise operators in model equations.
 
         Parameters
         ----------
         symbol : :class:`pybamm.expression_tree.symbol.Symbol` (or subclass) instance
             Symbol to discretise
-        y_slices : dict of {variable: slice}
-            The slices to assign to StateVectors when discretising a variable
-            (default None).
-        boundary_conditions : dict of {variable: boundary conditions}
-            Boundary conditions of the model
 
         Returns
         -------
@@ -204,49 +186,51 @@ class BaseDiscretisation(object):
 
         """
         if isinstance(symbol, pybamm.Gradient):
-            return self.gradient(symbol.children[0], y_slices, boundary_conditions)
+            child = symbol.children[0]
+            discretised_child = self.process_symbol(child)
+            return self._spatial_methods[symbol.domain[0]].gradient(
+                child, discretised_child, self._bcs
+            )
 
         elif isinstance(symbol, pybamm.Divergence):
-            return self.divergence(symbol.children[0], y_slices, boundary_conditions)
+            child = symbol.children[0]
+            discretised_child = self.process_symbol(child)
+            return self._spatial_methods[symbol.domain[0]].divergence(
+                child, discretised_child, self._bcs
+            )
 
         elif isinstance(symbol, pybamm.Broadcast):
             # Process child first
-            new_child = self.process_symbol(
-                symbol.children[0], y_slices, boundary_conditions
-            )
+            new_child = self.process_symbol(symbol.children[0])
             # Broadcast new_child to the domain specified by symbol.domain
             # Different discretisations may broadcast differently
-            return self.broadcast(new_child, symbol.domain)
+            if symbol.domain == []:
+                symbol = pybamm.NumpyBroadcast(symbol, symbol.domain, {})
+            else:
+                symbol = self._spatial_methods[symbol.domain[0]].broadcast(
+                    new_child, symbol.domain
+                )
+            return symbol
 
         elif isinstance(symbol, pybamm.BinaryOperator):
-            return self.process_binary_operators(symbol, y_slices, boundary_conditions)
+            return self.process_binary_operators(symbol)
 
         elif isinstance(symbol, pybamm.UnaryOperator):
-            new_child = self.process_symbol(
-                symbol.children[0], y_slices, boundary_conditions
-            )
+            new_child = self.process_symbol(symbol.children[0])
             return symbol.__class__(new_child)
 
         elif isinstance(symbol, pybamm.Variable):
-            assert isinstance(y_slices, dict), ValueError(
-                """y_slices should be dict, not {}""".format(type(y_slices))
-            )
-            return pybamm.StateVector(y_slices[symbol.id])
+            return pybamm.StateVector(self._y_slices[symbol.id])
 
         elif isinstance(symbol, pybamm.Space):
-            symbol_mesh = self.mesh.combine_submeshes(*symbol.domain)
-            return pybamm.Vector(symbol_mesh.nodes)
+            return self._spatial_methods[symbol.domain[0]].spatial_variable(symbol)
 
         elif isinstance(symbol, pybamm.Concatenation):
-            new_children = [
-                self.process_symbol(child, y_slices, boundary_conditions)
-                for child in symbol.children
-            ]
+            new_children = [self.process_symbol(child) for child in symbol.children]
             new_symbol = pybamm.DomainConcatenation(new_children, self.mesh)
 
             if new_symbol.is_constant():
                 return pybamm.Vector(new_symbol.evaluate())
-
             return new_symbol
 
         else:
@@ -254,7 +238,7 @@ class BaseDiscretisation(object):
             new_symbol.parent = None
             return new_symbol
 
-    def process_binary_operators(self, bin_op, y_slices, boundary_conditions):
+    def process_binary_operators(self, bin_op):
         """Discretise binary operators in model equations.  Performs appropriate
         averaging of diffusivities if one of the children is a gradient operator, so
         that discretised sizes match up.  This is mainly an issue for the Finite Volume
@@ -265,11 +249,6 @@ class BaseDiscretisation(object):
         ----------
         bin_op : :class:`pybamm.BinaryOperator` (or subclass)
             Binary operator to discretise
-        y_slices : dict of {variable: slice}
-            The slices to assign to StateVectors when discretising a variable
-            (default None).
-        boundary_conditions : dict of {variable: boundary conditions}
-            Boundary conditions of the model
 
         Returns
         -------
@@ -279,8 +258,8 @@ class BaseDiscretisation(object):
         """
         # Pre-process children
         left, right = bin_op.children
-        new_left = self.process_symbol(left, y_slices, boundary_conditions)
-        new_right = self.process_symbol(right, y_slices, boundary_conditions)
+        new_left = self.process_symbol(left)
+        new_right = self.process_symbol(right)
         # Post-processing to make sure discretised dimensions match
         # If neither child has gradients, or both children have gradients
         # no need to do any averaging
@@ -294,102 +273,46 @@ class BaseDiscretisation(object):
             left.has_gradient_and_not_divergence()
             and not right.has_gradient_and_not_divergence()
         ):
-            new_right = self.compute_diffusivity(new_right)
+            new_right = self._spatial_methods[bin_op.domain[0]].compute_diffusivity(
+                new_right
+            )
         # If only right child has gradient, compute diffusivity for left child
         elif (
             right.has_gradient_and_not_divergence()
             and not left.has_gradient_and_not_divergence()
         ):
-            new_left = self.compute_diffusivity(new_left)
+            new_left = self._spatial_methods[bin_op.domain[0]].compute_diffusivity(
+                new_left
+            )
         # Return new binary operator with appropriate class
         return bin_op.__class__(new_left, new_right)
 
-    def compute_diffusivity(self, symbol):
-        """Compute diffusivity; default behaviour is identity operator"""
-        return symbol
-
-    def gradient(self, symbol, y_slices, boundary_conditions):
-        """How to discretise gradient operators.
-
-        Parameters
-        ----------
-        symbol : :class:`Symbol` (or subclass)
-            The symbol (typically a variable) of which to take the gradient
-        y_slices : slice
-            The slice to assign to StateVector when discretising a variable
-        boundary_conditions : dict
-            The boundary conditions of the model
-            ({symbol.id: {"left": left bc, "right": right bc}})
-
-        """
-        raise NotImplementedError
-
-    def divergence(self, symbol, y_slices, boundary_conditions):
-        """How to discretise divergence operators.
-
-        Parameters
-        ----------
-        symbol : :class:`Symbol` (or subclass)
-            The symbol (typically a variable) of which to take the divergence
-        y_slices : slice
-            The slice to assign to StateVector when discretising a variable
-        boundary_conditions : dict
-            The boundary conditions of the model
-
-        """
-        raise NotImplementedError
-
-    def broadcast(self, symbol, domain):
-        """
-        Broadcast symbol to a specified domain. To do this, calls
-        :class:`pybamm.NumpyBroadcast`
-
-        Parameters
-        ----------
-        symbol : :class:`pybamm.Symbol`
-            The symbol to be broadcasted
-        domain : iterable of string
-            The domain to broadcast to
-        """
-        # create broadcasted symbol
-        broadcasted_symbol = pybamm.NumpyBroadcast(symbol, domain, self.mesh)
-
-        # if the broadcasted symbol evaluates to a constant value, replace the
-        # symbol-Vector multiplication with a single array
-        if broadcasted_symbol.is_constant():
-            broadcasted_symbol = pybamm.Array(
-                broadcasted_symbol.evaluate(), domain=broadcasted_symbol.domain
-            )
-
-        return broadcasted_symbol
-
     def concatenate(self, *symbols):
         return pybamm.NumpyModelConcatenation(*symbols)
 
-    def _concatenate_init(self, var_eqn_dict, y_slices):
+    def _concatenate_init(self, var_eqn_dict):
         """
         Concatenate a dictionary of {variable: equation} initial conditions using
-        y_slices
+        self._y_slices
 
-        The keys/variables in `var_eqn_dict` must be the same as the ids in `y_slices`
+        The keys/variables in `var_eqn_dict` must be the same as the ids in
+        `self._y_slices`.
         The resultant concatenation is ordered according to the ordering of the slice
-        values in `y_slices`
+        values in `self._y_slices`
 
         Parameters
         ----------
         var_eqn_dict : dict
             Equations ({variable: equation} dict) to dicretise
-        y_slices : dict of {variable id: slice}
-            The slices to assign to StateVectors when discretising
 
-        Returns
+                Returns
         -------
         var_eqn_dict : dict
             Discretised right-hand side equations
 
         """
         ids = {v.id for v in var_eqn_dict.keys()}
-        if ids != y_slices.keys():
+        if ids != self._y_slices.keys():
             given_variable_names = [v.name for v in var_eqn_dict.keys()]
             raise pybamm.ModelError(
                 "Initial conditions are insufficient. Only "
@@ -397,7 +320,7 @@ class BaseDiscretisation(object):
             )
 
         equations = list(var_eqn_dict.values())
-        slices = [y_slices[var.id] for var in var_eqn_dict.keys()]
+        slices = [self._y_slices[var.id] for var in var_eqn_dict.keys()]
 
         # sort equations according to slices
         sorted_equations = [eq for _, eq in sorted(zip(slices, equations))]
diff --git a/pybamm/expression_tree/concatenations.py b/pybamm/expression_tree/concatenations.py
index 3cd0b3e6..8a4940ed 100644
--- a/pybamm/expression_tree/concatenations.py
+++ b/pybamm/expression_tree/concatenations.py
@@ -53,7 +53,7 @@ class NumpyModelConcatenation(pybamm.Symbol):
     Upon evaluation, equations are concatenated using numpy concatenation.
     Unlike :class:`pybamm.Concatenation`, this doesn't check domains, as its only use
     is to concatenate model equations (e.g. rhs equations or initial conditions, in
-    :class:`pybamm.BaseDiscretisation`), which might have common domains
+    :class:`pybamm.Discretisation`), which might have common domains
 
     **Extends**: :class:`pybamm.Symbol`
 
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index f8d46916..17ce41fd 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -89,7 +89,7 @@ class SpatialOperator(UnaryOperator):
 
     Derived classes will specify the particular operator
 
-    This type of node will be replaced by the :class:`BaseDiscretisation`
+    This type of node will be replaced by the :class:`Discretisation`
     class with a :class:`Matrix`
 
     **Extends:** :class:`UnaryOperator`
@@ -172,19 +172,19 @@ class NumpyBroadcast(Broadcast):
         child node
     domain : iterable of string
         the domain to broadcast the child to
-    mesh : mesh class
-        the mesh used for discretisation
+    npts : dict
+        contains the number of points on each domain
 
     **Extends:** :class:`SpatialOperator`
     """
 
-    def __init__(self, child, domain, mesh):
+    def __init__(self, child, domain, npts):
         super().__init__(child, domain, name="numpy broadcast")
         # determine broadcasting vector size (size 1 if the domain is empty)
         if domain == []:
             self.broadcasting_vector_size = 1
         else:
-            self.broadcasting_vector_size = sum([mesh[dom].npts for dom in domain])
+            self.broadcasting_vector_size = sum([npts[dom] for dom in domain])
         # 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/variable.py b/pybamm/expression_tree/variable.py
index 0e3bdf0e..89e57763 100644
--- a/pybamm/expression_tree/variable.py
+++ b/pybamm/expression_tree/variable.py
@@ -9,7 +9,7 @@ import pybamm
 class Variable(pybamm.Symbol):
     """A node in the expression tree represending a dependent variable
 
-    This node will be discretised by :class:`.BaseDiscretisation` and converted
+    This node will be discretised by :class:`.Discretisation` and converted
     to a :class:`.Vector` node.
 
     Parameters
diff --git a/pybamm/meshes/__init__.py b/pybamm/meshes/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/pybamm/meshes/meshes.py b/pybamm/meshes/meshes.py
index f31ab1fe..14392cdf 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/meshes/meshes.py
@@ -37,6 +37,7 @@ class Mesh(dict):
 
     def __init__(self, geometry, submesh_types, submesh_pts):
         super().__init__()
+        self.submesh_pts = submesh_pts
         for domain in geometry:
             submesh_type = submesh_types[domain]
             submesh_pt = submesh_pts[domain]
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 78184343..b09a5337 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -54,6 +54,7 @@ class BaseModel(object):
         )
 
         self.default_geometry = pybamm.Geometry1DMacro()
+
         self.default_parameter_values.process_geometry(self.default_geometry)
         # provide mesh properties
         submesh_pts = {
@@ -67,8 +68,15 @@ class BaseModel(object):
             "positive electrode": pybamm.Uniform1DSubMesh,
         }
 
+        self.default_spatial_methods = {
+            "negative electrode": pybamm.FiniteVolume,
+            "separator": pybamm.FiniteVolume,
+            "positive electrode": pybamm.FiniteVolume,
+        }
         self.mesh = pybamm.Mesh(self.default_geometry, submesh_types, submesh_pts)
-        self.default_discretisation = pybamm.FiniteVolumeDiscretisation(self.mesh)
+        self.default_discretisation = pybamm.Discretisation(
+            self.mesh, self.default_spatial_methods
+        )
         self.default_solver = pybamm.ScipySolver(method="RK45")
 
     def _set_dict(self, dict, name):
diff --git a/pybamm/spatial_methods/__init__.py b/pybamm/spatial_methods/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/pybamm/discretisations/finite_volume_discretisations.py b/pybamm/spatial_methods/finite_volume.py
similarity index 77%
rename from pybamm/discretisations/finite_volume_discretisations.py
rename to pybamm/spatial_methods/finite_volume.py
index 7f50883d..b172ad3d 100644
--- a/pybamm/discretisations/finite_volume_discretisations.py
+++ b/pybamm/spatial_methods/finite_volume.py
@@ -9,52 +9,66 @@ import numpy as np
 from scipy.sparse import spdiags
 
 
-class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
-    """Discretisation using Finite Volumes.
+class FiniteVolume(pybamm.SpatialMethod):
+    """
+    A class which implements the steps specific to the finite volume method during
+    discretisation.
 
     Parameters
     ----------
-    mesh : :class:`pybamm.BaseMesh` (or subclass)
-        The underlying mesh for discretisation
+    mesh : :class:`pybamm.Mesh` (or subclass)
+        Contains all the submeshes for discretisation
 
-    **Extends:** :class:`pybamm.BaseDiscretisation`
+    **Extends:"": :class:`pybamm.SpatialMethod`
     """
 
     def __init__(self, mesh):
-        super().__init__(mesh)
+        self._mesh = mesh
 
-    def compute_diffusivity(self, symbol):
+    def spatial_variable(self, symbol):
         """
-        Compute the diffusivity at cell edges, based on the diffusivity at cell nodes.
-        For now we just take the arithemtic mean, though it may be better to take the
-        harmonic mean based on [1].
-
-        [1] Recktenwald, Gerald. "The control-volume finite-difference approximation to
-        the diffusion equation." (2012).
+        Creates a discretised spatial variable compatible with
+        the FiniteVolume method.
 
         Parameters
-        ----------
-        symbol : :class:`pybamm.Symbol`
-            Symbol to be averaged. When evaluated, this symbol returns either a scalar
-            or an array of shape (n,), where n is the number of points in the mesh for
-            the symbol's domain (n = self.mesh[symbol.domain].npts)
+        -----------
+        symbol : :class:`pybamm.SpatialVariable`
+            The spatial variable to be discretised.
 
         Returns
         -------
-        :class:`pybamm.NodeToEdge`
-            Averaged symbol. When evaluated, this returns either a scalar or an array of
-            shape (n-1,) as appropriate.
+        :class:`pybamm.Vector`
+            Contains the discretised spatial variable
         """
 
-        def arithmetic_mean(array):
-            """Calculate the arithemetic mean of an array"""
-            return (array[1:] + array[:-1]) / 2
+        # for finite volume we use the cell centres
+        symbol_mesh = self._mesh.combine_submeshes(*symbol.domain)
+        return pybamm.Vector(symbol_mesh.nodes)
 
-        return pybamm.NodeToEdge(symbol, arithmetic_mean)
+    def broadcast(self, symbol, domain):
+        """
+        Broadcast symbol to a specified domain. To do this, calls
+        :class:`pybamm.NumpyBroadcast`
+
+        See :meth: `pybamm.SpatialMethod.broadcast`
+        """
+
+        # for finite volume we send variables to cells and so use number_of_cells
+        number_of_cells = {dom: submesh.npts for dom, submesh in self._mesh.items()}
+        broadcasted_symbol = pybamm.NumpyBroadcast(symbol, domain, number_of_cells)
+
+        # if the broadcasted symbol evaluates to a constant value, replace the
+        # symbol-Vector multiplication with a single array
+        if broadcasted_symbol.is_constant():
+            broadcasted_symbol = pybamm.Array(
+                broadcasted_symbol.evaluate(), domain=broadcasted_symbol.domain
+            )
+
+        return broadcasted_symbol
 
-    def gradient(self, symbol, y_slices, boundary_conditions):
+    def gradient(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the gradient operator.
-        See :meth:`pybamm.BaseDiscretisation.gradient`
+        See :meth:`pybamm.SpatialMethod.gradient`
         """
         # Check that boundary condition keys are hashes (ids)
         for key in boundary_conditions.keys():
@@ -62,7 +76,6 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
                 "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:
@@ -79,54 +92,6 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
         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.
@@ -143,7 +108,7 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
             The (sparse) finite volume gradient matrix for the domain
         """
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self._mesh.combine_submeshes(*domain)
 
         # Create matrix using submesh
         n = submesh.npts
@@ -155,33 +120,33 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
         matrix = spdiags(data, diags, n - 1, n)
         return pybamm.Matrix(matrix)
 
-    def divergence(self, symbol, y_slices, boundary_conditions):
+    def divergence(self, symbol, discretised_symbol, boundary_conditions):
         """Matrix-vector multiplication to implement the divergence operator.
-        See :meth:`pybamm.BaseDiscretisation.gradient`
+        See :meth:`pybamm.SpatialMethod.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:
             # for the particles there will be a "negative particle" "left" and "right"
             # and also a "positive particle" left and right.
             lbc = boundary_conditions[symbol.id]["left"]
             rbc = boundary_conditions[symbol.id]["right"]
-            discretised_symbol = self.concatenate(lbc, discretised_symbol, rbc)
+            discretised_symbol = pybamm.NumpyModelConcatenation(
+                lbc, discretised_symbol, rbc
+            )
 
         domain = symbol.domain
-        # check for
+        # check for spherical domains
         if ("negative particle" or "positive particle") in domain:
 
             # implement spherical operator
             divergence_matrix = self.divergence_matrix(domain)
 
-            submesh = self.mesh.combine_submeshes(*domain)
+            submesh = self._mesh[domain[0]]
             r = pybamm.Vector(submesh.nodes)
             r_edges = pybamm.Vector(submesh.edges)
 
@@ -210,7 +175,7 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
             The (sparse) finite volume divergence matrix for the domain
         """
         # Create appropriate submesh by combining submeshes in domain
-        submesh = self.mesh.combine_submeshes(*domain)
+        submesh = self._mesh.combine_submeshes(*domain)
 
         # Create matrix using submesh
         n = submesh.npts + 1
@@ -222,6 +187,89 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
         matrix = spdiags(data, diags, n - 1, n)
         return pybamm.Matrix(matrix)
 
+    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 pybamm.NumpyModelConcatenation(
+            left_ghost_cell, discretised_symbol, right_ghost_cell
+        )
+
+    #######################################################
+    # Can probably be moved outside of the spatial method
+    ######################################################
+
+    def compute_diffusivity(self, symbol):
+        """
+        Compute the diffusivity at cell edges, based on the diffusivity at cell nodes.
+        For now we just take the arithemtic mean, though it may be better to take the
+        harmonic mean based on [1].
+
+        [1] Recktenwald, Gerald. "The control-volume finite-difference approximation to
+        the diffusion equation." (2012).
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.Symbol`
+            Symbol to be averaged. When evaluated, this symbol returns either a scalar
+            or an array of shape (n,), where n is the number of points in the mesh for
+            the symbol's domain (n = self.mesh[symbol.domain].npts)
+
+        Returns
+        -------
+        :class:`pybamm.NodeToEdge`
+            Averaged symbol. When evaluated, this returns either a scalar or an array of
+            shape (n-1,) as appropriate.
+        """
+
+        def arithmetic_mean(array):
+            """Calculate the arithemetic mean of an array"""
+            return (array[1:] + array[:-1]) / 2
+
+        return pybamm.NodeToEdge(symbol, arithmetic_mean)
+
 
 class NodeToEdge(pybamm.SpatialOperator):
     """A node in the expression tree representing a unary operator that evaluates the
diff --git a/pybamm/spatial_methods/spatial_method.py b/pybamm/spatial_methods/spatial_method.py
new file mode 100644
index 00000000..1006ed32
--- /dev/null
+++ b/pybamm/spatial_methods/spatial_method.py
@@ -0,0 +1,113 @@
+#
+# A general spatial method class
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+
+
+class SpatialMethod:
+    """
+    A general spatial methods class.
+    All spatial methods will follow the general form of SpatialMethod in
+    that they contain a method for broadcasting variables onto a mesh,
+    a gradient operator, and a diverence operator.
+
+    Parameters
+    ----------
+    mesh : :class: `pybamm.Mesh` (or subclass)
+        Contains all the submeshes for discretisation
+    """
+
+    def __init__(self, mesh):
+        self.mesh = mesh
+
+    def spatial_variable(self, symbol):
+        """
+        Creates a discretised spatial variable compatible with
+        the FiniteVolume method.
+
+        Parameters
+        -----------
+        symbol : :class:`pybamm.SpatialVariable`
+            The spatial variable to be discretised.
+
+        Returns
+        -------
+        :class:`pybamm.Vector`
+            Contains the discretised spatial variable
+        """
+        raise NotImplementedError
+
+    def broadcast(self, symbol, domain):
+        """
+        Broadcast symbol to a specified domain.
+
+        Parameters
+        ----------
+        symbol : :class:`pybamm.Symbol`
+            The symbol to be broadcasted
+        domain : iterable of string
+            The domain to broadcast to
+
+        Returns
+        -------
+        broadcasted_symbol: class: `pybamm.Array`
+            The discretised symbol of the correct size for
+            the spatial method
+        """
+        raise NotImplementedError
+
+    def gradient(self, symbol, discretised_symbol, boundary_conditions):
+        """
+        Implements the gradient for a spatial method.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol that we will take the gradient of.
+        broadcasted_symbol: class: pybamm.Array
+            The discretised symbol of the correct size
+
+        boundary_conditions : dict
+            The boundary conditions of the model
+            ({symbol.id: {"left": left bc, "right": right bc}})
+
+        Returns
+        -------
+        :class: `pybamm.Array`
+            Contains the result of acting the discretised gradient on
+            the child discretised_symbol
+        """
+        raise NotImplementedError
+
+    def divergence(self, symbol, broadcasted_symbol, boundary_conditions):
+        """
+        Implements the divergence for a spatial method.
+
+        Parameters
+        ----------
+        symbol: :class:`pybamm.Symbol`
+            The symbol that we will take the gradient of.
+        broadcasted_symbol: class: pybamm.Array
+            The discretised symbol of the correct size
+
+        boundary_conditions : dict
+            The boundary conditions of the model
+            ({symbol.id: {"left": left bc, "right": right bc}})
+        Returns
+        -------
+        :class: `pybamm.Array`
+            Contains the result of acting the discretised divergence on
+            the child discretised_symbol
+        """
+        raise NotImplementedError
+
+    # We could possibly move the following outside of SpatialMethod
+    # depending on the requirements of the FiniteVolume
+
+    def compute_diffusivity(self):
+        """Compute the diffusivity at edges of cells.
+        Could interpret this as: find diffusivity as
+        off grid locations
+        """
+        raise NotImplementedError


Test patch:
diff --git a/tests/shared.py b/tests/shared.py
index 20fc12ec..6fffc90d 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -37,6 +37,12 @@ class TestDefaults1DMacro:
 
         self.mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.submesh_pts)
 
+        self.spatial_methods = {
+            "negative electrode": SpatialMethodForTesting,
+            "separator": SpatialMethodForTesting,
+            "positive electrode": SpatialMethodForTesting,
+        }
+
 
 class TestDefaults1DParticle:
     def __init__(self, n):
@@ -52,25 +58,48 @@ class TestDefaults1DParticle:
 
         self.mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.submesh_pts)
 
+        self.spatial_methods = {"negative particle": pybamm.FiniteVolume}
 
-class DiscretisationForTesting(pybamm.BaseDiscretisation):
+
+class SpatialMethodForTesting(pybamm.SpatialMethod):
     """Identity operators, no boundary conditions."""
 
     def __init__(self, mesh):
+        self._mesh = mesh
         super().__init__(mesh)
 
-    def gradient(self, symbol, y_slices, boundary_conditions):
-        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+    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)
+
+    def broadcast(self, symbol, domain):
+        # for finite volume we send variables to cells and so use number_of_cells
+        number_of_cells = {dom: submesh.npts for dom, submesh in self._mesh.items()}
+        broadcasted_symbol = pybamm.NumpyBroadcast(symbol, domain, number_of_cells)
+
+        # if the broadcasted symbol evaluates to a constant value, replace the
+        # symbol-Vector multiplication with a single array
+        if broadcasted_symbol.is_constant():
+            broadcasted_symbol = pybamm.Array(
+                broadcasted_symbol.evaluate(), domain=broadcasted_symbol.domain
+            )
+
+        return broadcasted_symbol
+
+    def gradient(self, symbol, discretised_symbol, boundary_conditions):
         n = 0
         for domain in symbol.domain:
             n += self.mesh[domain].npts
         gradient_matrix = pybamm.Matrix(np.eye(n))
         return gradient_matrix * discretised_symbol
 
-    def divergence(self, symbol, y_slices, boundary_conditions):
-        discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
+    def divergence(self, symbol, discretised_symbol, boundary_conditions):
         n = 0
         for domain in symbol.domain:
             n += self.mesh[domain].npts
         divergence_matrix = pybamm.Matrix(np.eye(n))
         return divergence_matrix * discretised_symbol
+
+    def compute_diffusivity(self, symbol):
+        return symbol
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index 8f0379df..3eda96e3 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -13,7 +13,7 @@ class TestDiscretise(unittest.TestCase):
         a = pybamm.Variable("a")
         b = pybamm.Variable("b")
         c = pybamm.Variable("c")
-        y_slices = {c.id: slice(0, 1), a.id: slice(2, 3), b.id: slice(3, 4)}
+
         initial_conditions = {
             c: pybamm.Scalar(1),
             a: pybamm.Scalar(2),
@@ -22,9 +22,10 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
 
-        result = disc._concatenate_init(initial_conditions, y_slices)
+        disc._y_slices = {c.id: slice(0, 1), a.id: slice(2, 3), b.id: slice(3, 4)}
+        result = disc._concatenate_init(initial_conditions)
 
         self.assertIsInstance(result, pybamm.NumpyModelConcatenation)
         self.assertEqual(result.children[0].evaluate(), 1)
@@ -33,53 +34,57 @@ class TestDiscretise(unittest.TestCase):
 
         initial_conditions = {a: pybamm.Scalar(2), b: pybamm.Scalar(3)}
         with self.assertRaises(pybamm.ModelError):
-            result = disc._concatenate_init(initial_conditions, y_slices)
+            result = disc._concatenate_init(initial_conditions)
 
     def test_discretise_slicing(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        spatial_methods = {}
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         c = pybamm.Variable("c", domain=whole_cell)
-        variables = [c]
-        y_slices = disc.get_variable_slices(variables)
-        self.assertEqual(y_slices, {c.id: slice(0, 100)})
+        disc._variables = [c]
+        disc.set_variable_slices()
+
+        self.assertEqual(disc._y_slices, {c.id: slice(0, 100)})
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         c_true = combined_submesh.nodes ** 2
         y = c_true
-        np.testing.assert_array_equal(y[y_slices[c.id]], c_true)
+        np.testing.assert_array_equal(y[disc._y_slices[c.id]], c_true)
 
         # Several variables
         d = pybamm.Variable("d", domain=whole_cell)
         jn = pybamm.Variable("jn", domain=["negative electrode"])
-        variables = [c, d, jn]
-        y_slices = disc.get_variable_slices(variables)
+        disc._variables = [c, d, jn]
+        disc.set_variable_slices()
+
         self.assertEqual(
-            y_slices,
+            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
         y = np.concatenate([c_true, d_true, jn_true])
-        np.testing.assert_array_equal(y[y_slices[c.id]], c_true)
-        np.testing.assert_array_equal(y[y_slices[d.id]], d_true)
-        np.testing.assert_array_equal(y[y_slices[jn.id]], 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)
+        np.testing.assert_array_equal(y[disc._y_slices[jn.id]], jn_true)
 
     def test_process_symbol_base(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        spatial_methods = {}
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
 
         # variable
         var = pybamm.Variable("var")
-        y_slices = {var.id: slice(53)}
-        var_disc = disc.process_symbol(var, y_slices)
+        disc._y_slices = {var.id: slice(53)}
+        var_disc = disc.process_symbol(var)
         self.assertIsInstance(var_disc, pybamm.StateVector)
-        self.assertEqual(var_disc._y_slice, y_slices[var.id])
+        self.assertEqual(var_disc._y_slice, disc._y_slices[var.id])
         # scalar
         scal = pybamm.Scalar(5)
         scal_disc = disc.process_symbol(scal)
@@ -98,20 +103,20 @@ class TestDiscretise(unittest.TestCase):
 
         # binary operator
         bin = var + scal
-        bin_disc = disc.process_symbol(bin, y_slices)
+        bin_disc = disc.process_symbol(bin)
         self.assertIsInstance(bin_disc, pybamm.Addition)
         self.assertIsInstance(bin_disc.children[0], pybamm.StateVector)
         self.assertIsInstance(bin_disc.children[1], pybamm.Scalar)
 
         bin2 = scal + var
-        bin2_disc = disc.process_symbol(bin2, y_slices)
+        bin2_disc = disc.process_symbol(bin2)
         self.assertIsInstance(bin2_disc, pybamm.Addition)
         self.assertIsInstance(bin2_disc.children[0], pybamm.Scalar)
         self.assertIsInstance(bin2_disc.children[1], pybamm.StateVector)
 
         # non-spatial unary operator
         un1 = -var
-        un1_disc = disc.process_symbol(un1, y_slices)
+        un1_disc = disc.process_symbol(un1)
         self.assertIsInstance(un1_disc, pybamm.Negate)
         self.assertIsInstance(un1_disc.children[0], pybamm.StateVector)
 
@@ -131,10 +136,10 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
 
-        y_slices = {var1.id: slice(53), var2.id: slice(53, 59)}
-        exp_disc = disc.process_symbol(expression, y_slices)
+        disc._y_slices = {var1.id: slice(53), var2.id: slice(53, 59)}
+        exp_disc = disc.process_symbol(expression)
         self.assertIsInstance(exp_disc, pybamm.Division)
         # left side
         self.assertIsInstance(exp_disc.children[0], pybamm.Multiplication)
@@ -147,7 +152,8 @@ class TestDiscretise(unittest.TestCase):
             isinstance(exp_disc.children[0].children[1].children[1], pybamm.StateVector)
         )
         self.assertEqual(
-            exp_disc.children[0].children[1].children[1].y_slice, y_slices[var2.id]
+            exp_disc.children[0].children[1].children[1].y_slice,
+            disc._y_slices[var2.id],
         )
         # right side
         self.assertIsInstance(exp_disc.children[1], pybamm.Addition)
@@ -158,7 +164,8 @@ class TestDiscretise(unittest.TestCase):
             isinstance(exp_disc.children[1].children[0].children[0], pybamm.StateVector)
         )
         self.assertEqual(
-            exp_disc.children[1].children[0].children[0].y_slice, y_slices[var1.id]
+            exp_disc.children[1].children[0].children[0].y_slice,
+            disc._y_slices[var1.id],
         )
         self.assertTrue(
             isinstance(exp_disc.children[1].children[0].children[1], pybamm.Scalar)
@@ -167,17 +174,17 @@ class TestDiscretise(unittest.TestCase):
 
     def test_discretise_spatial_operator(self):
         # create discretisation
-        defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
-        mesh = disc.mesh
-
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=whole_cell)
-        y_slices = disc.get_variable_slices([var])
+        defaults = shared.TestDefaults1DMacro()
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        mesh = disc.mesh
+        disc._variables = [var]
+        disc.set_variable_slices()
 
         # Simple expressions
         for eqn in [pybamm.grad(var), pybamm.div(var)]:
-            eqn_disc = disc.process_symbol(eqn, y_slices, {})
+            eqn_disc = disc.process_symbol(eqn)
 
             self.assertIsInstance(eqn_disc, pybamm.Multiplication)
             self.assertIsInstance(eqn_disc.children[0], pybamm.Matrix)
@@ -185,7 +192,7 @@ class TestDiscretise(unittest.TestCase):
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
             y = combined_submesh.nodes ** 2
-            var_disc = disc.process_symbol(var, y_slices)
+            var_disc = disc.process_symbol(var)
             # grad and var are identity operators here (for testing purposes)
             np.testing.assert_array_equal(
                 eqn_disc.evaluate(None, y), var_disc.evaluate(None, y)
@@ -193,7 +200,7 @@ class TestDiscretise(unittest.TestCase):
 
         # More complex expressions
         for eqn in [var * pybamm.grad(var), var * pybamm.div(var)]:
-            eqn_disc = disc.process_symbol(eqn, y_slices, {})
+            eqn_disc = disc.process_symbol(eqn)
 
             self.assertIsInstance(eqn_disc, pybamm.Multiplication)
             self.assertIsInstance(eqn_disc.children[0], pybamm.StateVector)
@@ -202,7 +209,7 @@ class TestDiscretise(unittest.TestCase):
             self.assertIsInstance(eqn_disc.children[1].children[1], pybamm.StateVector)
 
             y = combined_submesh.nodes ** 2
-            var_disc = disc.process_symbol(var, y_slices)
+            var_disc = disc.process_symbol(var)
             # grad and var are identity operators here (for testing purposes)
             np.testing.assert_array_equal(
                 eqn_disc.evaluate(None, y), var_disc.evaluate(None, y) ** 2
@@ -210,12 +217,13 @@ class TestDiscretise(unittest.TestCase):
 
     def test_core_NotImplementedErrors(self):
         # create discretisation
-        disc = pybamm.BaseDiscretisation(None)
+        pybamm.Discretisation(None, {})
 
-        with self.assertRaises(NotImplementedError):
-            disc.gradient(None, None, {})
-        with self.assertRaises(NotImplementedError):
-            disc.divergence(None, None, {})
+        # TODO: implement and equivelent test in spatial_methods
+        # with self.assertRaises(NotImplementedError):
+        #     disc.gradient(None, None, {})
+        # with self.assertRaises(NotImplementedError):
+        #     disc.divergence(None, None, {})
 
     def test_process_dict(self):
         # one equation
@@ -225,19 +233,24 @@ class TestDiscretise(unittest.TestCase):
         rhs = {c: pybamm.div(N)}
         initial_conditions = {c: pybamm.Scalar(3)}
         variables = {"c_squared": c ** 2}
-        boundary_conditions = {N: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}}
+        boundary_conditions = {
+            N.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
+        }
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         y = combined_submesh.nodes ** 2
-        y_slices = disc.get_variable_slices(rhs.keys())
+        disc._variables = list(rhs.keys())
+        disc._bcs = boundary_conditions
+
+        disc.set_variable_slices()
         # rhs - grad and div are identity operators here
-        processed_rhs = disc.process_dict(rhs, y_slices, boundary_conditions)
+        processed_rhs = disc.process_dict(rhs)
         np.testing.assert_array_equal(y, processed_rhs[c].evaluate(None, y))
         # initial conditions
         y0 = disc.process_dict(initial_conditions)
@@ -245,7 +258,7 @@ class TestDiscretise(unittest.TestCase):
             y0[c].evaluate(0, None), 3 * np.ones_like(combined_submesh.nodes)
         )
         # vars
-        processed_vars = disc.process_dict(variables, y_slices, boundary_conditions)
+        processed_vars = disc.process_dict(variables)
         np.testing.assert_array_equal(
             y ** 2, processed_vars["c_squared"].evaluate(None, y)
         )
@@ -259,14 +272,16 @@ class TestDiscretise(unittest.TestCase):
         y = np.concatenate(
             [combined_submesh.nodes ** 2, mesh["negative electrode"].nodes ** 4]
         )
-        y_slices = disc.get_variable_slices(rhs.keys())
+
+        disc._variables = list(rhs.keys())
+        disc.set_variable_slices()
         # rhs
-        processed_rhs = disc.process_dict(rhs, y_slices, boundary_conditions)
+        processed_rhs = disc.process_dict(rhs)
         np.testing.assert_array_equal(
-            y[y_slices[c.id]], processed_rhs[c].evaluate(None, y)
+            y[disc._y_slices[c.id]], processed_rhs[c].evaluate(None, y)
         )
         np.testing.assert_array_equal(
-            y[y_slices[T.id]], processed_rhs[T].evaluate(None, y)
+            y[disc._y_slices[T.id]], processed_rhs[T].evaluate(None, y)
         )
         # initial conditions
         y0 = disc.process_dict(initial_conditions)
@@ -292,7 +307,7 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -371,7 +386,7 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         mesh = disc.mesh
 
         disc.process_model(model)
@@ -424,15 +439,20 @@ class TestDiscretise(unittest.TestCase):
     def test_broadcast(self):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
 
+        a = pybamm.Scalar(7)
+        vec = pybamm.Vector(np.linspace(0, 1))
+        var = pybamm.Variable("var")
+
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
+
         # scalar
-        a = pybamm.Scalar(7)
-        broad = disc.broadcast(a, whole_cell)
+        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)
@@ -440,8 +460,7 @@ class TestDiscretise(unittest.TestCase):
         self.assertEqual(broad.domain, whole_cell)
 
         # vector
-        vec = pybamm.Vector(np.linspace(0, 1))
-        broad = disc.broadcast(vec, ["separator"])
+        broad = disc._spatial_methods[whole_cell[0]].broadcast(vec, ["separator"])
         self.assertIsInstance(broad, pybamm.Array)
         np.testing.assert_array_equal(
             broad.evaluate(),
@@ -450,10 +469,9 @@ class TestDiscretise(unittest.TestCase):
         self.assertEqual(broad.domain, ["separator"])
 
         # process Broadcast symbol
-        var = pybamm.Variable("var")
-        y_slices = {var.id: slice(53)}
+        disc._y_slices = {var.id: slice(53)}
         broad1 = pybamm.Broadcast(var, ["negative electrode"])
-        broad1_disc = disc.process_symbol(broad1, y_slices)
+        broad1_disc = disc.process_symbol(broad1)
         self.assertIsInstance(broad1_disc, pybamm.NumpyBroadcast)
         self.assertIsInstance(broad1_disc.children[0], pybamm.StateVector)
 
@@ -470,7 +488,7 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
 
         conc = disc.concatenate(a, b, c)
         self.assertIsInstance(conc, pybamm.NumpyModelConcatenation)
@@ -482,14 +500,14 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         mesh = disc.mesh
 
-        var = pybamm.Variable("var", domain=whole_cell)
-        y_slices = disc.get_variable_slices([var])
+        disc._variables = [pybamm.Variable("var", domain=whole_cell)]
+        disc.set_variable_slices()
 
         eqn = pybamm.Concatenation(a, b)
-        eqn_disc = disc.process_symbol(eqn, y_slices, {})
+        eqn_disc = disc.process_symbol(eqn)
         self.assertIsInstance(eqn_disc, pybamm.Vector)
         expected_vector = np.concatenate(
             [
@@ -500,19 +518,22 @@ class TestDiscretise(unittest.TestCase):
         np.testing.assert_allclose(eqn_disc.evaluate(), expected_vector)
 
     def test_discretise_space(self):
+        # variables
+        x1 = pybamm.Space(["negative electrode"])
+        x2 = pybamm.Space(["negative electrode", "separator"])
+        x3 = 3 * pybamm.Space(["negative electrode"])
+
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
 
         # space
-        x1 = pybamm.Space(["negative electrode"])
         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
         )
 
-        x2 = pybamm.Space(["negative electrode", "separator"])
         x2_disc = disc.process_symbol(x2)
         self.assertIsInstance(x2_disc, pybamm.Vector)
         np.testing.assert_array_equal(
@@ -520,7 +541,6 @@ class TestDiscretise(unittest.TestCase):
             disc.mesh.combine_submeshes("negative electrode", "separator").nodes,
         )
 
-        x3 = 3 * pybamm.Space(["negative electrode"])
         x3_disc = disc.process_symbol(x3)
         self.assertIsInstance(x3_disc.children[1], pybamm.Vector)
         np.testing.assert_array_equal(
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index 1512db99..21cca07d 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -80,7 +80,7 @@ class TestConcatenations(unittest.TestCase):
     def test_numpy_domain_concatenation(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         mesh = disc.mesh
 
         a_dom = ["negative electrode"]
diff --git a/tests/test_expression_tree/test_unary_operators.py b/tests/test_expression_tree/test_unary_operators.py
index 7d3b44cd..304ba568 100644
--- a/tests/test_expression_tree/test_unary_operators.py
+++ b/tests/test_expression_tree/test_unary_operators.py
@@ -100,14 +100,15 @@ class TestUnaryOperators(unittest.TestCase):
     def test_numpy_broadcast(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         mesh = disc.mesh
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         combined_submeshes = mesh.combine_submeshes(*whole_cell)
         # scalar
         a = pybamm.Scalar(7)
-        broad = pybamm.NumpyBroadcast(a, whole_cell, mesh)
+        npts = {dom: submesh.npts for dom, submesh in mesh.items()}
+        broad = pybamm.NumpyBroadcast(a, whole_cell, npts)
         np.testing.assert_array_equal(
             broad.evaluate(), 7 * np.ones_like(combined_submeshes.nodes)
         )
@@ -115,7 +116,7 @@ class TestUnaryOperators(unittest.TestCase):
 
         # vector
         vec = pybamm.Vector(np.linspace(0, 1))
-        broad = pybamm.NumpyBroadcast(vec, whole_cell, mesh)
+        broad = pybamm.NumpyBroadcast(vec, whole_cell, npts)
         np.testing.assert_array_equal(
             broad.evaluate(),
             np.linspace(0, 1)[:, np.newaxis] * np.ones_like(combined_submeshes.nodes),
@@ -125,7 +126,7 @@ class TestUnaryOperators(unittest.TestCase):
 
         # state vector
         state_vec = pybamm.StateVector(slice(1, 2))
-        broad = pybamm.NumpyBroadcast(state_vec, whole_cell, mesh)
+        broad = pybamm.NumpyBroadcast(state_vec, whole_cell, npts)
         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
@@ -133,7 +134,7 @@ class TestUnaryOperators(unittest.TestCase):
 
         # state vector - bad input
         state_vec = pybamm.StateVector(slice(1, 5))
-        broad = pybamm.NumpyBroadcast(state_vec, whole_cell, mesh)
+        broad = pybamm.NumpyBroadcast(state_vec, whole_cell, npts)
         y = np.vstack([np.linspace(0, 1), np.linspace(0, 2)]).T
         with self.assertRaises(AssertionError):
             broad.evaluate(y=y)
diff --git a/tests/test_meshes/__init__.py b/tests/test_meshes/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/test_discretisations/test_meshes.py b/tests/test_meshes/test_meshes.py
similarity index 100%
rename from tests/test_discretisations/test_meshes.py
rename to tests/test_meshes/test_meshes.py
diff --git a/tests/test_meshes/test_submesh.py b/tests/test_meshes/test_submesh.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/test_meshes/test_uniform_submesh.py b/tests/test_meshes/test_uniform_submesh.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/test_models/test_interface.py b/tests/test_models/test_interface.py
index b83e2574..2b36824d 100644
--- a/tests/test_models/test_interface.py
+++ b/tests/test_models/test_interface.py
@@ -53,8 +53,8 @@ class TestHomogeneousReaction(unittest.TestCase):
         param = pybamm.ParameterValues(
             "input/parameters/lithium-ion/parameters/LCO.csv"
         )
-        mesh = shared.TestDefaults1DMacro().mesh
-        disc = pybamm.BaseDiscretisation(mesh)
+        defaults = shared.TestDefaults1DMacro()
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
 
         rxn = pybamm.interface.homogeneous_reaction()
 
@@ -162,11 +162,13 @@ class TestButlerVolmerLeadAcid(unittest.TestCase):
         param_bv_p = param.process_symbol(bv_p)
 
         # discretise
-        mesh = shared.TestDefaults1DMacro().mesh
-        disc = pybamm.BaseDiscretisation(mesh)
-        y_slices = disc.get_variable_slices([self.cn, self.cp, self.phin, self.phip])
-        processed_bv_n = disc.process_symbol(param_bv_n, y_slices)
-        processed_bv_p = disc.process_symbol(param_bv_p, y_slices)
+        defaults = shared.TestDefaults1DMacro()
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        mesh = disc.mesh
+        disc._variables = [self.cn, self.cp, self.phin, self.phip]
+        disc.set_variable_slices()
+        processed_bv_n = disc.process_symbol(param_bv_n)
+        processed_bv_p = disc.process_symbol(param_bv_p)
 
         submesh = np.concatenate(
             [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
@@ -194,10 +196,13 @@ class TestButlerVolmerLeadAcid(unittest.TestCase):
         param_bv_whole = param.process_symbol(bv_whole)
 
         # discretise
-        mesh = shared.TestDefaults1DMacro().mesh
-        disc = pybamm.BaseDiscretisation(mesh)
-        y_slices = disc.get_variable_slices([self.cn, self.cp, self.phin, self.phip])
-        processed_bv_whole = disc.process_symbol(param_bv_whole, y_slices)
+        defaults = shared.TestDefaults1DMacro()
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        mesh = disc.mesh
+
+        disc._variables = [self.cn, self.cp, self.phin, self.phip]
+        disc.set_variable_slices()
+        processed_bv_whole = disc.process_symbol(param_bv_whole)
 
         # test
         submesh = np.concatenate(
@@ -273,11 +278,14 @@ class TestExchangeCurrentDensity(unittest.TestCase):
         param_j0p = param.process_symbol(j0p)
 
         # discretise
-        mesh = shared.TestDefaults1DMacro().mesh
-        disc = pybamm.BaseDiscretisation(mesh)
-        y_slices = disc.get_variable_slices([self.cn, self.cp])
-        processed_j0n = disc.process_symbol(param_j0n, y_slices)
-        processed_j0p = disc.process_symbol(param_j0p, y_slices)
+        defaults = shared.TestDefaults1DMacro()
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
+        mesh = disc.mesh
+
+        disc._variables = [self.cn, self.cp]
+        disc.set_variable_slices()
+        processed_j0n = disc.process_symbol(param_j0n)
+        processed_j0p = disc.process_symbol(param_j0p)
 
         submesh = np.concatenate(
             [mesh["negative electrode"].nodes, mesh["positive electrode"].nodes]
diff --git a/tests/test_parameters/test_standard_parameters_lead_acid.py b/tests/test_parameters/test_standard_parameters_lead_acid.py
index c680b557..39363eb1 100644
--- a/tests/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/test_parameters/test_standard_parameters_lead_acid.py
@@ -60,8 +60,8 @@ class TestStandardParametersLeadAcid(unittest.TestCase):
         parameter_values = pybamm.ParameterValues(
             "input/parameters/lead-acid/default.csv", {"current scale": 1}
         )
-        mesh = shared.TestDefaults1DMacro().mesh
-        disc = pybamm.BaseDiscretisation(mesh)
+        defaults = shared.TestDefaults1DMacro()
+        disc = pybamm.Discretisation(defaults.mesh, defaults.spatial_methods)
         processed_s = disc.process_symbol(parameter_values.process_symbol(s))
 
         # test output
diff --git a/tests/test_solvers/test_scipy_solver.py b/tests/test_solvers/test_scipy_solver.py
index 599f3c84..469fe512 100644
--- a/tests/test_solvers/test_scipy_solver.py
+++ b/tests/test_solvers/test_scipy_solver.py
@@ -46,7 +46,12 @@ class TestScipySolver(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        spatial_methods = {
+            "negative electrode": pybamm.FiniteVolume,
+            "separator": pybamm.FiniteVolume,
+            "positive electrode": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         disc.process_model(model)
         # Solve
         solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
diff --git a/tests/test_spatial_methods/__init__.py b/tests/test_spatial_methods/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/tests/test_discretisations/test_finite_volume_discretisations.py b/tests/test_spatial_methods/test_finite_volume.py
similarity index 76%
rename from tests/test_discretisations/test_finite_volume_discretisations.py
rename to tests/test_spatial_methods/test_finite_volume.py
index 4e9cade2..7c7a2e71 100644
--- a/tests/test_discretisations/test_finite_volume_discretisations.py
+++ b/tests/test_spatial_methods/test_finite_volume.py
@@ -8,7 +8,7 @@ import numpy as np
 import unittest
 
 
-class TestFiniteVolumeDiscretisation(unittest.TestCase):
+class TestFiniteVolume(unittest.TestCase):
     def test_node_to_edge(self):
         a = pybamm.Symbol("a")
 
@@ -38,14 +38,20 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        spatial_methods = {
+            "negative electrode": pybamm.FiniteVolume,
+            "separator": pybamm.FiniteVolume,
+            "positive electrode": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         # Discretise some equations where averaging is needed
         var = pybamm.Variable("var", domain=whole_cell)
-        y_slices = disc.get_variable_slices([var])
+        disc._variables = [var]
+        disc.set_variable_slices()
         y_test = np.ones_like(combined_submesh.nodes)
         for eqn in [
             var * pybamm.grad(var),
@@ -57,7 +63,7 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             pybamm.grad(var) * pybamm.grad(var),
             (pybamm.grad(var) + 2) * pybamm.grad(var) ** 2,
         ]:
-            eqn_disc = disc.process_symbol(eqn, y_slices, {})
+            eqn_disc = disc.process_symbol(eqn)
             # Check that the equation can be evaluated
             eqn_disc.evaluate(None, y_test)
 
@@ -80,11 +86,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
                 -2 * pybamm.div(var * pybamm.grad(var) + 2 * pybamm.grad(var)),
             ],
         ):
-            eqn_disc = disc.process_symbol(
-                eqn,
-                y_slices,
-                {flux.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}},
-            )
+            disc._bcs = {flux.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}}
+
+            eqn_disc = disc.process_symbol(eqn)
             # Check that the equation can be evaluated
             eqn_disc.evaluate(None, y_test)
 
@@ -93,17 +97,25 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        spatial_methods = {
+            "negative electrode": pybamm.FiniteVolume,
+            "separator": pybamm.FiniteVolume,
+            "positive electrode": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         # Add ghost nodes
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         var = pybamm.Variable("var", domain=whole_cell)
-        y_slices = disc.get_variable_slices([var])
-        discretised_symbol = pybamm.StateVector(y_slices[var.id])
+        disc._variables = [var]
+        disc.set_variable_slices()
+        discretised_symbol = pybamm.StateVector(disc._y_slices[var.id])
         lbc = pybamm.Scalar(0)
         rbc = pybamm.Scalar(3)
-        symbol_plus_ghost = disc.add_ghost_nodes(discretised_symbol, lbc, rbc)
+        symbol_plus_ghost = pybamm.FiniteVolume(defaults.mesh).add_ghost_nodes(
+            discretised_symbol, lbc, rbc
+        )
 
         # Test
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -136,7 +148,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        spatial_methods = {
+            "negative electrode": pybamm.FiniteVolume,
+            "separator": pybamm.FiniteVolume,
+            "positive electrode": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -150,8 +167,11 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         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)
+        disc._bcs = boundary_conditions
+
+        disc._variables = [var]
+        disc.set_variable_slices()
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
         np.testing.assert_array_equal(
@@ -167,12 +187,14 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             var.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}
         }
 
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+        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)
         )
 
-        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        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)
         )
@@ -183,7 +205,8 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         """
         # create discretisation
         defaults = shared.TestDefaults1DParticle(10)
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        spatial_methods = {"negative particle": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes("negative particle")
@@ -198,8 +221,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         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)
+
+        disc._bcs = boundary_conditions
+
+        disc._variables = [var]
+        disc.set_variable_slices()
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
         np.testing.assert_array_equal(
@@ -210,8 +237,10 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         boundary_conditions = {
             var.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}
         }
+        disc._bcs = boundary_conditions
+
         y_linear = combined_submesh.nodes
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+        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)
         )
@@ -224,8 +253,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         boundary_conditions = {
             var.id: {"left": pybamm.Scalar(6), "right": pybamm.Scalar(6)}
         }
+        disc._bcs = boundary_conditions
 
-        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        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)
         )
@@ -236,7 +266,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        spatial_methods = {
+            "negative electrode": pybamm.FiniteVolume,
+            "separator": pybamm.FiniteVolume,
+            "positive electrode": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -244,8 +279,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         # grad
         var = pybamm.Variable("var", domain=whole_cell)
         grad_eqn = pybamm.grad(var)
-        y_slices = disc.get_variable_slices([var])
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, {})
+        disc._variables = [var]
+        disc.set_variable_slices()
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
         np.testing.assert_array_equal(
@@ -259,7 +295,8 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         boundary_conditions = {
             N.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
         }
-        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        disc._bcs = boundary_conditions
+        div_eqn_disc = disc.process_symbol(div_eqn)
 
         # Linear y should have laplacian zero
         linear_y = combined_submesh.nodes
@@ -276,7 +313,8 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DParticle(10)
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        spatial_methods = {"negative particle": pybamm.FiniteVolume}
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes("negative particle")
@@ -284,8 +322,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         # grad
         var = pybamm.Variable("var", domain="negative particle")
         grad_eqn = pybamm.grad(var)
-        y_slices = disc.get_variable_slices([var])
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, {})
+        disc._variables = [var]
+        disc.set_variable_slices()
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         constant_y = np.ones_like(combined_submesh.nodes)
         np.testing.assert_array_equal(
@@ -305,7 +344,8 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         boundary_conditions = {
             N.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
         }
-        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        disc._bcs = boundary_conditions
+        div_eqn_disc = disc.process_symbol(div_eqn)
 
         linear_y = combined_submesh.nodes
         const = 6 * np.ones(combined_submesh.npts)
@@ -320,7 +360,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         """
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        spatial_methods = {
+            "negative electrode": pybamm.FiniteVolume,
+            "separator": pybamm.FiniteVolume,
+            "positive electrode": pybamm.FiniteVolume,
+        }
+        disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
         mesh = disc.mesh
 
         mesh.add_ghost_meshes()
@@ -332,9 +377,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         boundary_conditions = {
             var.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
         }
-        y_slices = disc.get_variable_slices([var])
+        disc._bcs = boundary_conditions
+
+        disc._variables = [var]
+        disc.set_variable_slices()
 
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+        grad_eqn_disc = disc.process_symbol(grad_eqn)
 
         combined_submesh = mesh.combine_submeshes("negative electrode", "separator")
         constant_y = np.ones_like(combined_submesh.nodes)
@@ -353,13 +401,14 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
                 "right": pybamm.Scalar(combined_submesh.edges[-1]),
             }
         }
+        disc._bcs = boundary_conditions
 
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, 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)
         )
 
-        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
+        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)
         )
@@ -376,7 +425,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             n = 3 * round(n / 3)
             # create discretisation
             defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            spatial_methods = {
+                "negative electrode": pybamm.FiniteVolume,
+                "separator": pybamm.FiniteVolume,
+                "positive electrode": pybamm.FiniteVolume,
+            }
+            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
             mesh = disc.mesh
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -385,8 +439,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             grad_exact = np.cos(combined_submesh.edges[1:-1])
 
             # Discretise and evaluate
-            y_slices = disc.get_variable_slices([var])
-            grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, {})
+            disc._variables = [var]
+            disc.set_variable_slices()
+            grad_eqn_disc = disc.process_symbol(grad_eqn)
             grad_approx = grad_eqn_disc.evaluate(None, y)
 
             # Calculate errors
@@ -416,8 +471,14 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         def get_l2_error(n):
             # create discretisation
             defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            spatial_methods = {
+                "negative electrode": pybamm.FiniteVolume,
+                "separator": pybamm.FiniteVolume,
+                "positive electrode": pybamm.FiniteVolume,
+            }
+            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
             mesh = disc.mesh
+            disc._bcs = boundary_conditions
 
             # Set up discretisation
             whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -431,8 +492,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             grad_exact = np.cos(combined_submesh.edges)
 
             # Discretise and evaluate
-            y_slices = disc.get_variable_slices([var])
-            grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
+            disc._variables = [var]
+            disc.set_variable_slices()
+            grad_eqn_disc = disc.process_symbol(grad_eqn)
             grad_approx = grad_eqn_disc.evaluate(None, y)
 
             # Calculate errors
@@ -462,7 +524,13 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
 
             # create discretisation
             defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            spatial_methods = {
+                "negative electrode": pybamm.FiniteVolume,
+                "separator": pybamm.FiniteVolume,
+                "positive electrode": pybamm.FiniteVolume,
+            }
+            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
+            disc._bcs = boundary_conditions
             mesh = disc.mesh
 
             whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -473,8 +541,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             div_exact_internal = -np.sin(combined_submesh.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)
+            disc._variables = [var]
+            disc.set_variable_slices()
+            div_eqn_disc = disc.process_symbol(div_eqn)
             div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
 
             # Calculate errors
@@ -506,8 +575,14 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             whole_cell = ["negative electrode", "separator", "positive electrode"]
             # create discretisation
             defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            spatial_methods = {
+                "negative electrode": pybamm.FiniteVolume,
+                "separator": pybamm.FiniteVolume,
+                "positive electrode": pybamm.FiniteVolume,
+            }
+            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
             mesh = disc.mesh
+            disc._bcs = boundary_conditions
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
 
@@ -516,8 +591,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             div_exact_internal = -np.sin(combined_submesh.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)
+            disc._variables = [var]
+            disc.set_variable_slices()
+            div_eqn_disc = disc.process_symbol(div_eqn)
             div_approx_internal = div_eqn_disc.evaluate(None, y)[1:-1]
 
             # Calculate errors
@@ -548,7 +624,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         def get_l2_error(n):
             defaults = shared.TestDefaults1DParticle(n)
 
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            spatial_methods = {"negative particle": pybamm.FiniteVolume}
+            disc = pybamm.Discretisation(defaults.mesh, spatial_methods)
+            disc._bcs = boundary_conditions
             mesh = disc.mesh["negative particle"]
             r = mesh.nodes
 
@@ -558,8 +636,9 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             exact_internal = exact[1:-1]
 
             # discretise and evaluate
-            y_slices = disc.get_variable_slices([c])
-            eqn_disc = disc.process_symbol(eqn, y_slices, boundary_conditions)
+            disc._variables = [c]
+            disc.set_variable_slices()
+            eqn_disc = disc.process_symbol(eqn)
             approx_internal = eqn_disc.evaluate(None, y)[1:-1]
 
             # error