Repository: tinosulzer/PyBaMM

Issue:
add a `Function` node to expression tree
for `Parameters` that vary with time, user puts in the name of the function (functions defined in <name>.py in the same directory as the parameter file). Then the parameter class replaces `Parameter` with `Function` of that name

Discussion:

Patch:
diff --git a/docs/index.rst b/docs/index.rst
index 216607ff..e999b810 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -23,8 +23,6 @@ Contents
 .. toctree::
 
     source/expression_tree/index
-    source/geometry/index
-    source/meshes/index
     source/models/index
     source/parameters/index
     source/discretisations/index
diff --git a/docs/source/geometry/geometry.rst b/docs/source/geometry/geometry.rst
deleted file mode 100644
index d99e3948..00000000
--- a/docs/source/geometry/geometry.rst
+++ /dev/null
@@ -1,14 +0,0 @@
-Geometry
-==========
-
-.. autoclass:: pybamm.Geometry
-    :members:
-
-.. autoclass:: pybamm.Geometry1DMacro
-    :members:
-
-.. autoclass:: pybamm.Geometry1DMicro
-    :members:
-    
-.. autoclass:: pybamm.Geometry3DMacro
-    :members:
diff --git a/docs/source/geometry/index.rst b/docs/source/geometry/index.rst
deleted file mode 100644
index d2dc988d..00000000
--- a/docs/source/geometry/index.rst
+++ /dev/null
@@ -1,7 +0,0 @@
-Geometry
-========
-
-.. toctree::
-
-  geometry
- 
diff --git a/docs/source/meshes/index.rst b/docs/source/meshes/index.rst
deleted file mode 100644
index 7a518d02..00000000
--- a/docs/source/meshes/index.rst
+++ /dev/null
@@ -1,7 +0,0 @@
-Meshes
-======
-
-.. toctree::
-  
-  meshes
-  submeshes
diff --git a/docs/source/meshes/meshes.rst b/docs/source/meshes/meshes.rst
deleted file mode 100644
index f6788bb7..00000000
--- a/docs/source/meshes/meshes.rst
+++ /dev/null
@@ -1,5 +0,0 @@
-Meshes
-==========
-
-.. autoclass:: pybamm.Mesh
-    :members:
diff --git a/docs/source/meshes/submeshes.rst b/docs/source/meshes/submeshes.rst
deleted file mode 100644
index ea9d6b44..00000000
--- a/docs/source/meshes/submeshes.rst
+++ /dev/null
@@ -1,8 +0,0 @@
-Sub Meshes
-==========
-
-.. autoclass:: pybamm.SubMesh1D
-    :members:
-
-.. autoclass:: pybamm.Uniform1DSubMesh
-    :members:
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index 7bf5b0a3..89880d97 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -45,6 +45,7 @@ ABSOLUTE_PATH = os.path.join(os.path.split(script_path)[0], "..")
 #
 from .util import Timer
 from .util import profile
+from .util import load_function
 
 #
 # Classes for the Expression Tree
@@ -70,6 +71,7 @@ from .expression_tree.unary_operators import (
     UnaryOperator,
     Negate,
     AbsoluteValue,
+    Function,
     SpatialOperator,
     Gradient,
     Divergence,
@@ -78,6 +80,7 @@ from .expression_tree.unary_operators import (
     grad,
     div,
 )
+from .expression_tree.function_parameter import FunctionParameter
 from .expression_tree.scalar import Scalar
 from .expression_tree.variable import Variable
 from .expression_tree.independent_variable import IndependentVariable, Time, Space
@@ -125,9 +128,8 @@ from .discretisations.finite_volume_discretisations import (
     FiniteVolumeDiscretisation,
     NodeToEdge,
 )
-from .meshes.meshes import KNOWN_DOMAINS
-from .meshes.meshes import Mesh
-from .meshes.submeshes import SubMesh1D, Uniform1DSubMesh
+from .discretisations.meshes import KNOWN_DOMAINS
+from .discretisations.meshes import Mesh, SubMesh1D, Uniform1DSubMesh
 
 #
 # Simulation class
diff --git a/pybamm/discretisations/base_discretisation.py b/pybamm/discretisations/base_discretisation.py
index 067cfdb7..fb66ed60 100644
--- a/pybamm/discretisations/base_discretisation.py
+++ b/pybamm/discretisations/base_discretisation.py
@@ -5,7 +5,6 @@ from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
 
-import copy
 import numpy as np
 
 
@@ -23,14 +22,25 @@ class BaseDiscretisation(object):
         the types of submeshes being employed on each of the subdomains
     """
 
-    def __init__(self, mesh):
-        self._mesh = mesh
+    def __init__(self, mesh_type, submesh_pts, submesh_types):
+        self.mesh_type = mesh_type
+        self.submesh_pts = submesh_pts
+        self.submesh_types = submesh_types
+
+        self._mesh = {}
 
     @property
     def mesh(self):
         return self._mesh
 
-    def process_model(self, model):
+    @mesh.setter
+    def mesh(self, mesh_in):
+        self._mesh = mesh_in
+
+    def mesh_geometry(self, geometry):
+        self._mesh = self.mesh_type(geometry, self.submesh_types, self.submesh_pts)
+
+    def process_model(self, model, geometry):
         """Discretise a model.
         Currently inplace, could be changed to return a new model.
 
@@ -42,6 +52,9 @@ class BaseDiscretisation(object):
 
         """
 
+        # Set the mesh
+        self._mesh = self.mesh_type(geometry, self.submesh_types, self.submesh_pts)
+
         # Set the y split for variables
         variables = list(model.rhs.keys()) + list(model.algebraic.keys())
         y_slices = self.get_variable_slices(variables)
@@ -249,10 +262,22 @@ class BaseDiscretisation(object):
 
             return new_symbol
 
+        elif isinstance(symbol, pybamm.Scalar):
+            return pybamm.Scalar(symbol.value, domain=symbol.domain)
+
+        elif isinstance(symbol, pybamm.Array):
+            return symbol.__class__(symbol.entries, domain=symbol.domain)
+
         else:
-            new_symbol = copy.deepcopy(symbol)
-            new_symbol.parent = None
-            return new_symbol
+            raise NotImplementedError
+            # hack to copy the symbol but without a parent
+            # (building tree from bottom up)
+            # simply setting new_symbol.parent = None, after copying, raises a TreeError
+            # parent = symbol.parent
+            # symbol.parent = None
+            # new_symbol = copy.copy(symbol)
+            # symbol.parent = parent
+            # return new_symbol
 
     def process_binary_operators(self, bin_op, y_slices, boundary_conditions):
         """Discretise binary operators in model equations.  Performs appropriate
diff --git a/pybamm/discretisations/finite_volume_discretisations.py b/pybamm/discretisations/finite_volume_discretisations.py
index 7f50883d..f02fc70c 100644
--- a/pybamm/discretisations/finite_volume_discretisations.py
+++ b/pybamm/discretisations/finite_volume_discretisations.py
@@ -20,8 +20,8 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
     **Extends:** :class:`pybamm.BaseDiscretisation`
     """
 
-    def __init__(self, mesh):
-        super().__init__(mesh)
+    def __init__(self, mesh_type, submesh_pts, submesh_types):
+        super().__init__(mesh_type, submesh_pts, submesh_types)
 
     def compute_diffusivity(self, symbol):
         """
@@ -74,8 +74,6 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
                 + domain
                 + [domain[-1] + "_right ghost cell"]
             )
-
-        # note in 1D spherical grad and normal grad are the same
         gradient_matrix = self.gradient_matrix(domain)
         return gradient_matrix * discretised_symbol
 
@@ -168,31 +166,11 @@ class FiniteVolumeDiscretisation(pybamm.BaseDiscretisation):
         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)
-
-        domain = symbol.domain
-        # check for
-        if ("negative particle" or "positive particle") in domain:
-
-            # implement spherical operator
-            divergence_matrix = self.divergence_matrix(domain)
-
-            submesh = self.mesh.combine_submeshes(*domain)
-            r = pybamm.Vector(submesh.nodes)
-            r_edges = pybamm.Vector(submesh.edges)
-
-            out = (1 / (r ** 2)) * (
-                divergence_matrix * ((r_edges ** 2) * discretised_symbol)
-            )
-
-        else:
-            divergence_matrix = self.divergence_matrix(domain)
-            out = divergence_matrix * discretised_symbol
-        return out
+        divergence_matrix = self.divergence_matrix(symbol.domain)
+        return divergence_matrix * discretised_symbol
 
     def divergence_matrix(self, domain):
         """
diff --git a/pybamm/meshes/meshes.py b/pybamm/discretisations/meshes.py
similarity index 67%
rename from pybamm/meshes/meshes.py
rename to pybamm/discretisations/meshes.py
index f31ab1fe..49a51303 100644
--- a/pybamm/meshes/meshes.py
+++ b/pybamm/discretisations/meshes.py
@@ -7,14 +7,7 @@ from __future__ import print_function, unicode_literals
 import numpy as np
 import pybamm
 
-KNOWN_DOMAINS = [
-    "negative electrode",
-    "separator",
-    "positive electrode",
-    "test",
-    "negative particle",
-    "positive particle",
-]
+KNOWN_DOMAINS = ["negative electrode", "separator", "positive electrode", "test"]
 
 
 class Mesh(dict):
@@ -74,7 +67,7 @@ class Mesh(dict):
                 [self[submeshnames[0]].edges]
                 + [self[submeshname].edges[1:] for submeshname in submeshnames[1:]]
             )
-            return pybamm.SubMesh1D(combined_submesh_edges)
+            return SubMesh1D(combined_submesh_edges)
         else:
             raise pybamm.DomainError("submesh edges are not aligned")
 
@@ -93,8 +86,47 @@ class Mesh(dict):
 
             # left ghost cell: two edges, one node, to the left of existing submesh
             lgs_edges = np.array([2 * edges[0] - edges[1], edges[0]])
-            self[submeshname + "_left ghost cell"] = pybamm.SubMesh1D(lgs_edges)
+            self[submeshname + "_left ghost cell"] = SubMesh1D(lgs_edges)
 
             # right ghost cell: two edges, one node, to the right of existing submesh
             rgs_edges = np.array([edges[-1], 2 * edges[-1] - edges[-2]])
-            self[submeshname + "_right ghost cell"] = pybamm.SubMesh1D(rgs_edges)
+            self[submeshname + "_right ghost cell"] = SubMesh1D(rgs_edges)
+
+
+class SubMesh1D:
+    """ Submesh class.
+        Contains the position of the nodes and the number of mesh points.
+
+        Parameters
+        ----------
+        domain : dict
+            A dictionary that contains the limits of the spatial variables
+        npts : dict
+            A dictionary that contains the number of points to be used on each
+            spatial variable
+        """
+
+    def __init__(self, edges):
+        self.edges = edges
+        self.nodes = (self.edges[1:] + self.edges[:-1]) / 2
+        self.d_edges = np.diff(self.edges)
+        self.d_nodes = np.diff(self.nodes)
+        self.npts = self.nodes.size
+
+
+class Uniform1DSubMesh(SubMesh1D):
+    """
+    A class to generate a uniform submesh on a 1D domain
+    """
+
+    def __init__(self, lims, npts):
+
+        # currently accept lims and npts as dicts. This may get changed at a future
+        # date depending on the form of mesh we desire for 2D/3D
+
+        spatial_lims = list(lims.values())[0]
+        npts = list(npts.values())[0]
+
+        edges = np.linspace(spatial_lims["min"], spatial_lims["max"], npts + 1)
+
+        super().__init__(edges)
diff --git a/pybamm/expression_tree/function_parameter.py b/pybamm/expression_tree/function_parameter.py
new file mode 100644
index 00000000..7cc3a1cd
--- /dev/null
+++ b/pybamm/expression_tree/function_parameter.py
@@ -0,0 +1,27 @@
+#
+# Parameter class
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+import pybamm
+
+
+class FunctionParameter(pybamm.UnaryOperator):
+    """A node in the expression tree representing a function parameter
+
+    This node will be replaced by a :class:`.Function` node
+
+    Parameters
+    ----------
+
+    name : str
+        name of the node
+    child : :class:`Symbol`
+        child node
+    domain : iterable of str, optional
+        list of domains the parameter is valid over, defaults to empty list
+
+    """
+
+    def __init__(self, name, child):
+        super().__init__(name, child)
diff --git a/pybamm/expression_tree/symbol.py b/pybamm/expression_tree/symbol.py
index 78d4ef01..d88662bb 100644
--- a/pybamm/expression_tree/symbol.py
+++ b/pybamm/expression_tree/symbol.py
@@ -36,6 +36,7 @@ class Symbol(anytree.NodeMixin):
             # copy child before adding
             # this also adds copy.copy(child) to self.children
             copy.copy(child).parent = self
+
         self.domain = domain
 
     @property
@@ -100,19 +101,6 @@ class Symbol(anytree.NodeMixin):
             + tuple(self.domain)
         )
 
-    @property
-    def orphans(self):
-        """
-        Returning deepcopies of the children, with parents removed to avoid corrupting
-        the expression tree internal data
-        """
-        orp = []
-        for child in self.children:
-            new_child = copy.deepcopy(child)
-            new_child.parent = None
-            orp.append(new_child)
-        return tuple(orp)
-
     def render(self):
         """print out a visual representation of the tree (this node and its
         children)
diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py
index e23b8317..f8d46916 100644
--- a/pybamm/expression_tree/unary_operators.py
+++ b/pybamm/expression_tree/unary_operators.py
@@ -67,6 +67,22 @@ class AbsoluteValue(UnaryOperator):
         return np.abs(self.children[0].evaluate(t, y))
 
 
+class Function(UnaryOperator):
+    """A node in the expression tree representing an arbitrary function
+
+    **Extends:** :class:`UnaryOperator`
+    """
+
+    def __init__(self, func, child):
+        """ See :meth:`pybamm.UnaryOperator.__init__()`. """
+        super().__init__("function ({})".format(func.__name__), child)
+        self.func = func
+
+    def evaluate(self, t=None, y=None):
+        """ See :meth:`pybamm.Symbol.evaluate()`. """
+        return self.func(self.children[0].evaluate(t, y))
+
+
 class SpatialOperator(UnaryOperator):
     """A node in the expression tree representing a unary spatial operator
     (e.g. grad, div)
diff --git a/pybamm/geometry/geometry.py b/pybamm/geometry/geometry.py
index 8b47ee8c..62a3e2f8 100644
--- a/pybamm/geometry/geometry.py
+++ b/pybamm/geometry/geometry.py
@@ -25,14 +25,12 @@ class Geometry(dict):
 class Geometry1DMacro(Geometry):
     def __init__(self, custom_geometry={}):
         super().__init__()
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        x = pybamm.IndependentVariable("x", whole_cell)
         ln = pybamm.standard_parameters.ln
         ls = pybamm.standard_parameters.ls
 
-        self["negative electrode"] = {x: {"min": pybamm.Scalar(0), "max": ln}}
-        self["separator"] = {x: {"min": ln, "max": ln + ls}}
-        self["positive electrode"] = {x: {"min": ln + ls, "max": pybamm.Scalar(1)}}
+        self["negative electrode"] = {"x": {"min": pybamm.Scalar(0), "max": ln}}
+        self["separator"] = {"x": {"min": ln, "max": ln + ls}}
+        self["positive electrode"] = {"x": {"min": ln + ls, "max": pybamm.Scalar(1)}}
 
         # update with custom geometry if non empty
         self.update(custom_geometry)
@@ -42,13 +40,11 @@ class Geometry1DMicro(Geometry):
     def __init__(self, custom_geometry={}):
         super().__init__()
 
-        r = pybamm.IndependentVariable("r", ["negative particle", "positive particle"])
-
         self["negative particle"] = {
-            r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
+            "r": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
         }
         self["positive particle"] = {
-            r: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
+            "r": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
         }
         # update with custom geometry if non empty
         self.update(custom_geometry)
@@ -68,16 +64,11 @@ class Geometry3DMacro(Geometry1DMacro):
     def __init__(self, custom_geometry={}):
         super().__init__()
 
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        y = pybamm.IndependentVariable("y", whole_cell)
-        z = pybamm.IndependentVariable("z", whole_cell)
-
-        y_lim = {"min": pybamm.Scalar(0), "max": pybamm.standard_parameters.ly}
+        y = {"min": pybamm.Scalar(0), "max": pybamm.standard_parameters.ly}
 
-        z_lim = {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
+        z = {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
 
-        MACRO_DOMAINS = ["negative electrode", "separator", "positive electrode"]
-        for domain in MACRO_DOMAINS:
-            self[domain][y] = y_lim
-            self[domain][z] = z_lim
+        for domain in self["macro"]:
+            self["macro"][domain]["y"] = y
+            self["macro"][domain]["z"] = z
         self.update(custom_geometry)
diff --git a/pybamm/meshes/submeshes.py b/pybamm/meshes/submeshes.py
deleted file mode 100644
index 9914b305..00000000
--- a/pybamm/meshes/submeshes.py
+++ /dev/null
@@ -1,46 +0,0 @@
-#
-# Native PyBaMM Meshes
-#
-from __future__ import absolute_import, division
-from __future__ import print_function, unicode_literals
-
-import numpy as np
-
-
-class SubMesh1D:
-    """ Submesh class.
-        Contains the position of the nodes and the number of mesh points.
-
-        Parameters
-        ----------
-        domain : dict
-            A dictionary that contains the limits of the spatial variables
-        npts : dict
-            A dictionary that contains the number of points to be used on each
-            spatial variable
-        """
-
-    def __init__(self, edges):
-        self.edges = edges
-        self.nodes = (self.edges[1:] + self.edges[:-1]) / 2
-        self.d_edges = np.diff(self.edges)
-        self.d_nodes = np.diff(self.nodes)
-        self.npts = self.nodes.size
-
-
-class Uniform1DSubMesh(SubMesh1D):
-    """
-    A class to generate a uniform submesh on a 1D domain
-    """
-
-    def __init__(self, lims, npts):
-
-        # currently accept lims and npts as dicts. This may get changed at a future
-        # date depending on the form of mesh we desire for 2D/3D
-
-        spatial_lims = list(lims.values())[0]
-        npts = list(npts.values())[0]
-
-        edges = np.linspace(spatial_lims["min"], spatial_lims["max"], npts + 1)
-
-        super().__init__(edges)
diff --git a/pybamm/models/base_model.py b/pybamm/models/base_model.py
index 78184343..d79435a2 100644
--- a/pybamm/models/base_model.py
+++ b/pybamm/models/base_model.py
@@ -67,8 +67,11 @@ class BaseModel(object):
             "positive electrode": pybamm.Uniform1DSubMesh,
         }
 
-        self.mesh = pybamm.Mesh(self.default_geometry, submesh_types, submesh_pts)
-        self.default_discretisation = pybamm.FiniteVolumeDiscretisation(self.mesh)
+        mesh_type = pybamm.Mesh
+
+        self.default_discretisation = pybamm.FiniteVolumeDiscretisation(
+            mesh_type, submesh_pts, submesh_types
+        )
         self.default_solver = pybamm.ScipySolver(method="RK45")
 
     def _set_dict(self, dict, name):
diff --git a/pybamm/parameters/parameter_values.py b/pybamm/parameters/parameter_values.py
index eeb60a94..c2d199c0 100644
--- a/pybamm/parameters/parameter_values.py
+++ b/pybamm/parameters/parameter_values.py
@@ -146,6 +146,11 @@ class ParameterValues(dict):
             value = self.get_parameter_value(symbol)
             return pybamm.Scalar(value, domain=symbol.domain)
 
+        elif isinstance(symbol, pybamm.FunctionParameter):
+            new_child = self.process_symbol(symbol.children[0])
+            function_name = self.get_parameter_value(symbol)
+            return pybamm.Function(pybamm.load_function(function_name), new_child)
+
         elif isinstance(symbol, pybamm.BinaryOperator):
             left, right = symbol.children
             new_left = self.process_symbol(left)
@@ -168,6 +173,11 @@ class ParameterValues(dict):
             return pybamm.Concatenation(*new_children)
 
         else:
-            new_symbol = copy.deepcopy(symbol)
-            new_symbol.parent = None
+            # hack to copy the symbol but without a parent
+            # (building tree from bottom up)
+            # simply setting new_symbol.parent = None, after copying, raises a TreeError
+            parent = symbol.parent
+            symbol.parent = None
+            new_symbol = copy.copy(symbol)
+            symbol.parent = parent
             return new_symbol
diff --git a/pybamm/parameters/standard_parameters.py b/pybamm/parameters/standard_parameters.py
index d7117f46..efa65d97 100644
--- a/pybamm/parameters/standard_parameters.py
+++ b/pybamm/parameters/standard_parameters.py
@@ -82,10 +82,6 @@ Ls = Parameter("Ls")
 Lp = Parameter("Lp")
 Lx = Ln + Ls + Lp
 
-# 3D Geometry
-Ly = Parameter("Ly")
-Lz = Parameter("Lz")
-
 # Microscale Geometry
 R_n = Parameter("R_n")
 R_p = Parameter("R_p")
@@ -123,8 +119,6 @@ cp0_dimensional = Parameter("cp0")  # Initial li concentration in pos electrode
 ln = Ln / Lx
 ls = Ls / Lx
 lp = Lp / Lx
-ly = Ly / Lz
-lz = Lz / Lz
 
 # Microscale Geometry
 epsilon_n = Parameter("epsilon_n")  # Electrolyte volume fraction in neg electrode
diff --git a/pybamm/util.py b/pybamm/util.py
index a6f34f66..bd4fca3f 100644
--- a/pybamm/util.py
+++ b/pybamm/util.py
@@ -7,9 +7,12 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 
-import timeit
 import cProfile
+import importlib
+import os
 import pstats
+import sys
+import timeit
 
 
 class Timer(object):
@@ -78,3 +81,74 @@ def profile(code, sort="cumulative", num=30):
     stats = pstats.Stats()
     stats.sort_stats(sort)
     return stats
+
+
+def load_function(filename):
+    """
+    Load a python function from a file "function_name.py" called "function_name".
+    The filename might either be an absolute path, in which case that specific file will
+    be used, or the file will be searched for relative to PyBaMM root.
+
+    Arguments
+    ---------
+    filename : str
+        The name of the file containing the function of the same name.
+
+    Returns
+    -------
+    function
+        The python function loaded from the file.
+    """
+
+    if not filename.endswith('.py'):
+        raise ValueError("Expected filename.py, but got {}".format(filename))
+
+    # If it's an absolute path, find that exact file
+    if os.path.isabs(filename):
+        if not os.path.isfile(filename):
+            raise ValueError(
+                "{} is an absolute path, but the file is not found".format(filename))
+
+        valid_filename = filename
+
+    # Else, search in the whole PyBaMM directory for matches
+    else:
+        search_path = os.path.commonpath([pth for pth in sys.path if 'PyBaMM' in pth])
+        head, tail = os.path.split(filename)
+
+        matching_files = []
+
+        for root, _, files in os.walk(search_path):
+            for file in files:
+                if file == tail:
+                    full_path = os.path.join(root, file)
+                    if full_path.endswith(filename):
+                        matching_files.append(full_path)
+
+        if len(matching_files) == 0:
+            raise ValueError(
+                "{} cannot be found in the PyBaMM directory".format(filename))
+        elif len(matching_files) > 1:
+            raise ValueError(
+                "{} found multiple times in the PyBaMM directory".format(filename))
+
+        valid_filename = matching_files[0]
+
+    # Now: we have some /path/to/valid/filename.py
+    # Add "/path/to/vaid" to the python path, and load the module "filename".
+    # Then, check "filename" module contains "filename" function.  If it does, return
+    # that function object, or raise an exception
+
+    valid_path, valid_leaf = os.path.split(valid_filename)
+    sys.path.append(valid_path)
+
+    # Load the module, which must be the leaf of filename, minus the .py extension
+    valid_module = valid_leaf.replace('.py', '')
+    module_object = importlib.import_module(valid_module)
+
+    # Check that a function of the same name exists in the loaded module
+    if valid_module not in dir(module_object):
+        raise ValueError(
+            "No function {} found in module {}".format(valid_module, valid_module))
+
+    return getattr(module_object, valid_module)


Test patch:
diff --git a/tests/shared.py b/tests/shared.py
index 20fc12ec..6d1d86fc 100644
--- a/tests/shared.py
+++ b/tests/shared.py
@@ -7,7 +7,7 @@ import numpy as np
 
 
 class TestDefaults1DMacro:
-    def __init__(self, npts=0):
+    def __init__(self):
         self.param = pybamm.ParameterValues(
             base_parameters={"Ln": 0.3, "Ls": 0.3, "Lp": 0.3}
         )
@@ -20,44 +20,31 @@ class TestDefaults1DMacro:
             "separator": {"x": 25},
             "positive electrode": {"x": 35},
         }
-
         self.submesh_types = {
             "negative electrode": pybamm.Uniform1DSubMesh,
             "separator": pybamm.Uniform1DSubMesh,
             "positive electrode": pybamm.Uniform1DSubMesh,
         }
 
-        if npts != 0:
-            n = 3 * round(npts / 3)
-            self.submesh_pts = {
-                "negative electrode": {"x": n},
-                "separator": {"x": n},
-                "positive electrode": {"x": n},
-            }
-
-        self.mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.submesh_pts)
+        self.mesh_type = pybamm.Mesh
 
-
-class TestDefaults1DParticle:
-    def __init__(self, n):
-        self.geometry = {
-            "negative particle": {
-                "r": {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}
-            }
+    def set_equal_pts(self, npts):
+        n = 3 * round(npts / 3)
+        self.submesh_pts = {
+            "negative electrode": {"x": n},
+            "separator": {"x": n},
+            "positive electrode": {"x": n},
         }
-        self.param = pybamm.ParameterValues(base_parameters={})
-        self.param.process_geometry(self.geometry)
-        self.submesh_pts = {"negative particle": {"r": n}}
-        self.submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
-
-        self.mesh = pybamm.Mesh(self.geometry, self.submesh_types, self.submesh_pts)
 
 
 class DiscretisationForTesting(pybamm.BaseDiscretisation):
     """Identity operators, no boundary conditions."""
 
-    def __init__(self, mesh):
-        super().__init__(mesh)
+    def __init__(self, mesh_type, submesh_pts, submesh_types):
+        super().__init__(mesh_type, submesh_pts, submesh_types)
+
+    def set_mesh(self, geometry):
+        self._mesh = self.mesh_type(geometry, self.submesh_types, self.submesh_pts)
 
     def gradient(self, symbol, y_slices, boundary_conditions):
         discretised_symbol = self.process_symbol(symbol, y_slices, boundary_conditions)
diff --git a/tests/test_discretisations/test_base_discretisations.py b/tests/test_discretisations/test_base_discretisations.py
index 8f0379df..d24f0c94 100644
--- a/tests/test_discretisations/test_base_discretisations.py
+++ b/tests/test_discretisations/test_base_discretisations.py
@@ -22,7 +22,10 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
 
         result = disc._concatenate_init(initial_conditions, y_slices)
 
@@ -38,7 +41,10 @@ class TestDiscretise(unittest.TestCase):
     def test_discretise_slicing(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -72,7 +78,10 @@ class TestDiscretise(unittest.TestCase):
     def test_process_symbol_base(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
 
         # variable
         var = pybamm.Variable("var")
@@ -120,6 +129,10 @@ class TestDiscretise(unittest.TestCase):
         self.assertIsInstance(un2_disc, pybamm.AbsoluteValue)
         self.assertIsInstance(un2_disc.children[0], pybamm.Scalar)
 
+        # parameter should fail
+        with self.assertRaises(NotImplementedError):
+            disc.process_symbol(pybamm.Parameter("par"))
+
     def test_process_complex_expression(self):
         var1 = pybamm.Variable("var1")
         var2 = pybamm.Variable("var2")
@@ -131,7 +144,10 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
 
         y_slices = {var1.id: slice(53), var2.id: slice(53, 59)}
         exp_disc = disc.process_symbol(expression, y_slices)
@@ -168,7 +184,10 @@ class TestDiscretise(unittest.TestCase):
     def test_discretise_spatial_operator(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -210,7 +229,11 @@ class TestDiscretise(unittest.TestCase):
 
     def test_core_NotImplementedErrors(self):
         # create discretisation
-        disc = pybamm.BaseDiscretisation(None)
+        defaults = shared.TestDefaults1DMacro()
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
 
         with self.assertRaises(NotImplementedError):
             disc.gradient(None, None, {})
@@ -229,7 +252,10 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -292,11 +318,13 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.process_model(model, defaults.geometry)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
-        disc.process_model(model)
 
         y0 = model.concatenated_initial_conditions
         np.testing.assert_array_equal(y0, 3 * np.ones_like(combined_submesh.nodes))
@@ -324,7 +352,7 @@ class TestDiscretise(unittest.TestCase):
         }
         model.variables = {"ST": S * T}
 
-        disc.process_model(model)
+        disc.process_model(model, defaults.geometry)
         y0 = model.concatenated_initial_conditions
         np.testing.assert_array_equal(
             y0,
@@ -350,7 +378,7 @@ class TestDiscretise(unittest.TestCase):
         model.boundary_conditions = {}
         model.variables = {"ST": S * T}
         with self.assertRaises(pybamm.ModelError):
-            disc.process_model(model)
+            disc.process_model(model, defaults.geometry)
 
     def test_process_model_dae(self):
         # one rhs equation and one algebraic
@@ -371,10 +399,12 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.process_model(model, defaults.geometry)
         mesh = disc.mesh
 
-        disc.process_model(model)
         combined_submesh = mesh.combine_submeshes(*whole_cell)
 
         y0 = model.concatenated_initial_conditions
@@ -419,14 +449,17 @@ class TestDiscretise(unittest.TestCase):
         model.variables = {"c": c, "N": N, "d": d}
 
         with self.assertRaises(pybamm.ModelError):
-            disc.process_model(model)
+            disc.process_model(model, defaults.geometry)
 
     def test_broadcast(self):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -470,7 +503,10 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
 
         conc = disc.concatenate(a, b, c)
         self.assertIsInstance(conc, pybamm.NumpyModelConcatenation)
@@ -482,7 +518,10 @@ class TestDiscretise(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         var = pybamm.Variable("var", domain=whole_cell)
@@ -502,7 +541,10 @@ class TestDiscretise(unittest.TestCase):
     def test_discretise_space(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
 
         # space
         x1 = pybamm.Space(["negative electrode"])
diff --git a/tests/test_discretisations/test_finite_volume_discretisations.py b/tests/test_discretisations/test_finite_volume_discretisations.py
index 4e9cade2..ab5b4c56 100644
--- a/tests/test_discretisations/test_finite_volume_discretisations.py
+++ b/tests/test_discretisations/test_finite_volume_discretisations.py
@@ -38,7 +38,10 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        disc = pybamm.FiniteVolumeDiscretisation(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -93,7 +96,10 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        disc = pybamm.FiniteVolumeDiscretisation(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         # Add ghost nodes
@@ -136,7 +142,10 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        disc = pybamm.FiniteVolumeDiscretisation(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -177,66 +186,16 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             div_eqn_disc.evaluate(None, linear_y), np.zeros_like(combined_submesh.nodes)
         )
 
-    def test_spherical_grad_div_shapes_Dirichlet_bcs(self):
-        """
-        Test grad and div with Dirichlet boundary conditions (applied by grad on var)
-        """
-        # create discretisation
-        defaults = shared.TestDefaults1DParticle(10)
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
-        mesh = disc.mesh
-
-        combined_submesh = mesh.combine_submeshes("negative particle")
-
-        mesh.add_ghost_meshes()
-        disc.mesh.add_ghost_meshes()
-
-        # grad
-        # grad(r) == 1
-        var = pybamm.Variable("var", domain=["negative particle"])
-        grad_eqn = pybamm.grad(var)
-        boundary_conditions = {
-            var.id: {"left": pybamm.Scalar(1), "right": pybamm.Scalar(1)}
-        }
-        y_slices = disc.get_variable_slices([var])
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
-
-        constant_y = np.ones_like(combined_submesh.nodes)
-        np.testing.assert_array_equal(
-            grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges),
-        )
-
-        boundary_conditions = {
-            var.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(1)}
-        }
-        y_linear = combined_submesh.nodes
-        grad_eqn_disc = disc.process_symbol(grad_eqn, y_slices, boundary_conditions)
-        np.testing.assert_array_almost_equal(
-            grad_eqn_disc.evaluate(None, y_linear), np.ones_like(combined_submesh.edges)
-        )
-
-        # div: test on linear r^2
-        # div (grad r^2) = 6
-        const = 6 * np.ones(combined_submesh.npts)
-        N = pybamm.grad(var)
-        div_eqn = pybamm.div(N)
-        boundary_conditions = {
-            var.id: {"left": pybamm.Scalar(6), "right": pybamm.Scalar(6)}
-        }
-
-        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
-        np.testing.assert_array_almost_equal(
-            div_eqn_disc.evaluate(None, const), np.zeros_like(combined_submesh.nodes)
-        )
-
     def test_grad_div_shapes_Neumann_bcs(self):
         """Test grad and div with Neumann boundary conditions (applied by div on N)"""
         whole_cell = ["negative electrode", "separator", "positive electrode"]
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        disc = pybamm.FiniteVolumeDiscretisation(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -271,56 +230,16 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             div_eqn_disc.evaluate(None, linear_y), np.zeros_like(combined_submesh.nodes)
         )
 
-    def test_spherical_grad_div_shapes_Neumann_bcs(self):
-        """Test grad and div with Neumann boundary conditions (applied by div on N)"""
-
-        # create discretisation
-        defaults = shared.TestDefaults1DParticle(10)
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
-        mesh = disc.mesh
-
-        combined_submesh = mesh.combine_submeshes("negative particle")
-
-        # 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, {})
-
-        constant_y = np.ones_like(combined_submesh.nodes)
-        np.testing.assert_array_equal(
-            grad_eqn_disc.evaluate(None, constant_y),
-            np.zeros_like(combined_submesh.edges[1:-1]),
-        )
-
-        linear_y = combined_submesh.nodes
-        np.testing.assert_array_almost_equal(
-            grad_eqn_disc.evaluate(None, linear_y),
-            np.ones_like(combined_submesh.edges[1:-1]),
-        )
-        # div
-        # div ( grad(r^2) ) == 6 , N_left = N_right = 0
-        N = pybamm.grad(var)
-        div_eqn = pybamm.div(N)
-        boundary_conditions = {
-            N.id: {"left": pybamm.Scalar(0), "right": pybamm.Scalar(0)}
-        }
-        div_eqn_disc = disc.process_symbol(div_eqn, y_slices, boundary_conditions)
-
-        linear_y = combined_submesh.nodes
-        const = 6 * np.ones(combined_submesh.npts)
-
-        np.testing.assert_array_almost_equal(
-            div_eqn_disc.evaluate(None, const), np.zeros_like(combined_submesh.nodes)
-        )
-
     def test_grad_div_shapes_mixed_domain(self):
         """
         Test grad and div with Dirichlet boundary conditions (applied by grad on var)
         """
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        disc = pybamm.FiniteVolumeDiscretisation(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         mesh.add_ghost_meshes()
@@ -375,8 +294,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
             # Set up discretisation
             n = 3 * round(n / 3)
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            defaults = shared.TestDefaults1DMacro()
+            defaults.set_equal_pts(n)
+            disc = pybamm.FiniteVolumeDiscretisation(
+                defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+            )
+            disc.mesh_geometry(defaults.geometry)
             mesh = disc.mesh
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -415,8 +338,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         # Function for convergence testing
         def get_l2_error(n):
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            defaults = shared.TestDefaults1DMacro()
+            defaults.set_equal_pts(n)
+            disc = pybamm.FiniteVolumeDiscretisation(
+                defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+            )
+            disc.mesh_geometry(defaults.geometry)
             mesh = disc.mesh
 
             # Set up discretisation
@@ -461,8 +388,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         def get_l2_error(n):
 
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            defaults = shared.TestDefaults1DMacro()
+            defaults.set_equal_pts(n)
+            disc = pybamm.FiniteVolumeDiscretisation(
+                defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+            )
+            disc.mesh_geometry(defaults.geometry)
             mesh = disc.mesh
 
             whole_cell = ["negative electrode", "separator", "positive electrode"]
@@ -505,8 +436,12 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         def get_l2_error(n):
             whole_cell = ["negative electrode", "separator", "positive electrode"]
             # create discretisation
-            defaults = shared.TestDefaults1DMacro(n)
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+            defaults = shared.TestDefaults1DMacro()
+            defaults.set_equal_pts(n)
+            disc = pybamm.FiniteVolumeDiscretisation(
+                defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+            )
+            disc.mesh_geometry(defaults.geometry)
             mesh = disc.mesh
 
             combined_submesh = mesh.combine_submeshes(*whole_cell)
@@ -534,48 +469,6 @@ class TestFiniteVolumeDiscretisation(unittest.TestCase):
         rates = np.log2(errs[:-1] / errs[1:])
         np.testing.assert_array_less(1.49 * np.ones_like(rates), rates)
 
-    def test_spherical_operators(self):
-        # test div( grad( sin(r) )) == (2/r)*cos(r) - *sin(r)
-
-        domain = ["negative particle"]
-        c = pybamm.Variable("c", domain=domain)
-        N = pybamm.grad(c)
-        eqn = pybamm.div(N)
-        boundary_conditions = {
-            N.id: {"left": pybamm.Scalar(np.cos(0)), "right": pybamm.Scalar(np.cos(1))}
-        }
-
-        def get_l2_error(n):
-            defaults = shared.TestDefaults1DParticle(n)
-
-            disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
-            mesh = disc.mesh["negative particle"]
-            r = mesh.nodes
-
-            # exact solution
-            y = np.sin(r)
-            exact = (2 / r) * np.cos(r) - np.sin(r)
-            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)
-            approx_internal = eqn_disc.evaluate(None, y)[1:-1]
-
-            # error
-            error = np.linalg.norm(approx_internal - exact_internal) / np.linalg.norm(
-                exact_internal
-            )
-            return error
-
-        # Get errors
-        ns = 10 * (2 ** np.arange(2, 7))
-        errs = np.array([get_l2_error(int(n)) for n in ns])
-
-        # Get rates: expect h**1.5 convergence because of boundary conditions
-        rates = np.log2(errs[:-1] / errs[1:])
-        np.testing.assert_array_less(1.99 * np.ones_like(rates), rates)
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_expression_tree/test_concatenations.py b/tests/test_expression_tree/test_concatenations.py
index 6043a316..04225588 100644
--- a/tests/test_expression_tree/test_concatenations.py
+++ b/tests/test_expression_tree/test_concatenations.py
@@ -80,7 +80,10 @@ class TestConcatenations(unittest.TestCase):
     def test_numpy_domain_concatenation(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         a_dom = ["negative electrode"]
@@ -125,24 +128,6 @@ class TestConcatenations(unittest.TestCase):
             ),
         )
 
-    def test_concatenation_orphans(self):
-        a = pybamm.Variable("a")
-        b = pybamm.Variable("b")
-        c = pybamm.Variable("c")
-        conc = pybamm.Concatenation(a, b, c)
-        a_new, b_new, c_new = conc.orphans
-
-        # We should be able to manipulate the children without TreeErrors
-        self.assertIsInstance(2 * a_new, pybamm.Multiplication)
-        self.assertIsInstance(3 + b_new, pybamm.Addition)
-        self.assertIsInstance(4 - c_new, pybamm.Subtraction)
-
-        # ids should stay the same
-        self.assertEqual(a.id, a_new.id)
-        self.assertEqual(b.id, b_new.id)
-        self.assertEqual(c.id, c_new.id)
-        self.assertEqual(conc.id, pybamm.Concatenation(a_new, b_new, c_new).id)
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_expression_tree/test_symbol.py b/tests/test_expression_tree/test_symbol.py
index 85b145ed..98b593f3 100644
--- a/tests/test_expression_tree/test_symbol.py
+++ b/tests/test_expression_tree/test_symbol.py
@@ -231,17 +231,6 @@ class TestSymbol(unittest.TestCase):
         self.assertFalse(algebraic_eqn.has_gradient())
         self.assertFalse(algebraic_eqn.has_divergence())
 
-    def test_orphans(self):
-        a = pybamm.Symbol("a")
-        b = pybamm.Symbol("b")
-        sum = a + b
-
-        a_orp, b_orp = sum.orphans
-        self.assertIsNone(a_orp.parent)
-        self.assertIsNone(b_orp.parent)
-        self.assertEqual(a.id, a_orp.id)
-        self.assertEqual(b.id, b_orp.id)
-
 
 if __name__ == "__main__":
     print("Add -v for more debug output")
diff --git a/tests/test_expression_tree/test_unary_operators.py b/tests/test_expression_tree/test_unary_operators.py
index f3bb3929..f5638487 100644
--- a/tests/test_expression_tree/test_unary_operators.py
+++ b/tests/test_expression_tree/test_unary_operators.py
@@ -9,6 +9,10 @@ import unittest
 import numpy as np
 
 
+def test_function(arg):
+    return arg + arg
+
+
 class TestUnaryOperators(unittest.TestCase):
     def test_unary_operator(self):
         a = pybamm.Symbol("a", domain=["test"])
@@ -36,6 +40,26 @@ class TestUnaryOperators(unittest.TestCase):
         absb = pybamm.AbsoluteValue(b)
         self.assertEqual(absb.evaluate(), 4)
 
+    def test_function(self):
+        a = pybamm.Symbol("a")
+        funca = pybamm.Function(test_function, a)
+        self.assertEqual(funca.name, "function (test_function)")
+        self.assertEqual(funca.children[0].name, a.name)
+
+        b = pybamm.Scalar(1)
+        sina = pybamm.Function(np.sin, b)
+        self.assertEqual(sina.evaluate(), np.sin(1))
+        self.assertEqual(sina.name, "function ({})".format(np.sin.__name__))
+
+        c = pybamm.Vector(np.linspace(0, 1))
+        cosb = pybamm.Function(np.cos, c)
+        np.testing.assert_array_equal(cosb.evaluate(), np.cos(c.evaluate()))
+
+        var = pybamm.StateVector(slice(0, 100))
+        y = np.linspace(0, 1, 100)
+        logvar = pybamm.Function(np.log1p, var)
+        np.testing.assert_array_equal(logvar.evaluate(y=y), np.log1p(y))
+
     def test_gradient(self):
         a = pybamm.Symbol("a")
         grad = pybamm.Gradient(a)
@@ -76,7 +100,10 @@ class TestUnaryOperators(unittest.TestCase):
     def test_numpy_broadcast(self):
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.mesh_geometry(defaults.geometry)
         mesh = disc.mesh
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
diff --git a/tests/test_geometry/test_geometry.py b/tests/test_geometry/test_geometry.py
index cf7236ea..e8df45e0 100644
--- a/tests/test_geometry/test_geometry.py
+++ b/tests/test_geometry/test_geometry.py
@@ -2,78 +2,19 @@
 # Tests for the base model class
 #
 import pybamm
-import unittest
 
 
-class TestGeometry1DMacro(unittest.TestCase):
+class TestGeometry:
     def test_add_custom_geometry(self):
-        geometry = pybamm.Geometry1DMacro()
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        x = pybamm.IndependentVariable("x", whole_cell)
         custom_geometry = {
             "negative electrode": {
-                x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
-            }
-        }
-        geometry.update(custom_geometry)
-        self.assertEqual(
-            geometry["negative electrode"], custom_geometry["negative electrode"]
-        )
-
-    def test_geometry_keys(self):
-        geometry = pybamm.Geometry1DMacro()
-        for spatial_vars in geometry.values():
-            all(
-                self.assertIsInstance(spatial_var, pybamm.IndependentVariable)
-                for spatial_var in spatial_vars.keys()
-            )
-
-
-class TestGeometry1DMicro(unittest.TestCase):
-    def test_add_custom_geometry(self):
-        geometry = pybamm.Geometry1DMicro()
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        x = pybamm.IndependentVariable("x", whole_cell)
-        custom_geometry = {
-            "negative electrode": {
-                x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
-            }
-        }
-        geometry.update(custom_geometry)
-
-        self.assertEqual(
-            geometry["negative electrode"], custom_geometry["negative electrode"]
-        )
-
-    def test_geometry_keys(self):
-        geometry = pybamm.Geometry1DMicro()
-        for spatial_vars in geometry.values():
-            all(
-                self.assertIsInstance(spatial_var, pybamm.IndependentVariable)
-                for spatial_var in spatial_vars.keys()
-            )
-
-
-class TestGeometry3DMacro(unittest.TestCase):
-    def test_add_custom_geometry(self):
-        geometry = pybamm.Geometry3DMacro()
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-        x = pybamm.IndependentVariable("x", whole_cell)
-        custom_geometry = {
-            "negative electrode": {
-                x: {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
+                "x": {"min": pybamm.Scalar(1), "max": pybamm.Scalar(2)}
             }
         }
+        geometry = pybamm.Geometry1DMacro(custom_geometry)
 
-        geometry.update(custom_geometry)
         self.assertEqual(
             geometry["negative electrode"], custom_geometry["negative electrode"]
         )
 
-    def test_geometry_keys(self):
-        geometry = pybamm.Geometry3DMacro()
-        for spatial_vars in geometry.values():
-            all(
-                self.assertIsInstance(spatial_var, pybamm.IndependentVariable)
-                for spatial_var in spatial_vars.keys()
-            )
+    # process geometry test conducted within TestParameterValues
diff --git a/tests/test_models/standard_model_tests.py b/tests/test_models/standard_model_tests.py
index 2c303aa2..4236cd80 100644
--- a/tests/test_models/standard_model_tests.py
+++ b/tests/test_models/standard_model_tests.py
@@ -27,7 +27,7 @@ class StandardModelTest(object):
         # Overwrite discretisation if given
         if disc is not None:
             self.disc = disc
-        self.disc.process_model(self.model)
+        self.disc.process_model(self.model, self.geometry)
         # Model should still be well-posed after processing
         self.model.check_well_posedness()
 
diff --git a/tests/test_models/test_interface.py b/tests/test_models/test_interface.py
index 5c3a1db2..9d6dd7b0 100644
--- a/tests/test_models/test_interface.py
+++ b/tests/test_models/test_interface.py
@@ -4,7 +4,6 @@
 from __future__ import absolute_import, division
 from __future__ import print_function, unicode_literals
 import pybamm
-import tests.shared as shared
 
 import unittest
 
@@ -50,16 +49,32 @@ class TestHomogeneousReaction(unittest.TestCase):
         self.assertEqual(processed_rxn.children[2].evaluate() * lp.evaluate(), -1)
 
     def test_discretisation(self):
-        defaults = shared.TestDefaults1DMacro()
-        disc = pybamm.FiniteVolumeDiscretisation(defaults.mesh)
+        whole_cell = ["negative electrode", "separator", "positive electrode"]
+        param = pybamm.ParameterValues(
+            "input/parameters/lithium-ion/parameters/LCO.csv"
+        )
+        geometry = pybamm.Geometry1DMacro()
+        param.process_geometry(geometry)
 
+        mesh_type = pybamm.Mesh
+        submesh_pts = {
+            "negative electrode": {"x": 40},
+            "separator": {"x": 25},
+            "positive electrode": {"x": 35},
+        }
+        submesh_types = {
+            "negative electrode": pybamm.Uniform1DSubMesh,
+            "separator": pybamm.Uniform1DSubMesh,
+            "positive electrode": pybamm.Uniform1DSubMesh,
+        }
+
+        disc = pybamm.FiniteVolumeDiscretisation(mesh_type, submesh_pts, submesh_types)
+        disc._mesh = mesh_type(geometry, submesh_types, submesh_pts)
         rxn = pybamm.interface.homogeneous_reaction()
 
-        param_rxn = defaults.param.process_symbol(rxn)
+        param_rxn = param.process_symbol(rxn)
         processed_rxn = disc.process_symbol(param_rxn)
 
-        whole_cell = ["negative electrode", "separator", "positive electrode"]
-
         combined_submeshes = disc.mesh.combine_submeshes(*whole_cell)
         # processed_rxn should be a vector with the right shape
         self.assertIsInstance(processed_rxn, pybamm.Vector)
diff --git a/tests/test_models/test_lead_acid.py b/tests/test_models/test_lead_acid.py
index 9677f20d..2ccb0e60 100644
--- a/tests/test_models/test_lead_acid.py
+++ b/tests/test_models/test_lead_acid.py
@@ -23,7 +23,7 @@ class TestLeadAcidLOQS(unittest.TestCase):
         # process parameter values, discretise and solve
         model.default_parameter_values.process_model(model)
         disc = model.default_discretisation
-        disc.process_model(model)
+        disc.process_model(model, model.default_geometry)
         t_eval = np.linspace(0, 1, 100)
         solver = model.default_solver
         solver.solve(model, t_eval)
diff --git a/tests/test_models/test_simple_ode_model.py b/tests/test_models/test_simple_ode_model.py
index b6523c80..904a1dd9 100644
--- a/tests/test_models/test_simple_ode_model.py
+++ b/tests/test_models/test_simple_ode_model.py
@@ -23,8 +23,9 @@ class TestSimpleODEModel(unittest.TestCase):
 
         # discretise and solve
         disc = model.default_discretisation
+        disc.mesh_geometry(model.default_geometry)
         combined_submesh = disc.mesh.combine_submeshes(*whole_cell)
-        disc.process_model(model)
+        disc.process_model(model, model.default_geometry)
         t_eval = np.linspace(0, 1, 100)
         solver = model.default_solver
         solver.solve(model, t_eval)
diff --git a/tests/test_parameters/data/process_symbol_bad_function.py b/tests/test_parameters/data/process_symbol_bad_function.py
new file mode 100644
index 00000000..dec6e336
--- /dev/null
+++ b/tests/test_parameters/data/process_symbol_bad_function.py
@@ -0,0 +1,10 @@
+def no_valid_function():
+    """
+    This is a test function for load_function in util.py.
+
+    load_function expects a file some_name.py to contain a function called some_name.
+    This file has a function with a different name to the filename, and this is used in
+    test_util.py to verify this case is handled correctly.
+
+    """
+    return None
diff --git a/tests/test_parameters/data/process_symbol_test_function.py b/tests/test_parameters/data/process_symbol_test_function.py
new file mode 100644
index 00000000..5ed5288a
--- /dev/null
+++ b/tests/test_parameters/data/process_symbol_test_function.py
@@ -0,0 +1,11 @@
+def process_symbol_test_function(var):
+    """
+    This function tests process_symbol for Function nodes.
+
+    Parameters
+    ----------
+    var: numeric
+        Some variable
+
+   """
+    return 123 * var
diff --git a/tests/test_parameters/test_parameter_values.py b/tests/test_parameters/test_parameter_values.py
index 7fc95f2a..0ae3b722 100644
--- a/tests/test_parameters/test_parameter_values.py
+++ b/tests/test_parameters/test_parameter_values.py
@@ -118,6 +118,22 @@ class TestParameterValues(unittest.TestCase):
         self.assertIsInstance(processed_f, pybamm.Matrix)
         np.testing.assert_array_equal(processed_f.evaluate(), np.ones((5, 6)))
 
+    def test_process_function_parameter(self):
+        parameter_values = pybamm.ParameterValues(
+            {"a": 3, "func": 'process_symbol_test_function.py'})
+
+        # process parameter
+        a = pybamm.Parameter("a")
+        processed_a = parameter_values.process_symbol(a)
+        self.assertIsInstance(processed_a, pybamm.Scalar)
+        self.assertEqual(processed_a.value, 3)
+
+        # process function
+        func = pybamm.FunctionParameter("func", a)
+        processed_func = parameter_values.process_symbol(func)
+        self.assertIsInstance(processed_func, pybamm.Function)
+        self.assertEqual(processed_func.evaluate(), 369)
+
     def test_process_complex_expression(self):
         var1 = pybamm.Variable("var1")
         var2 = pybamm.Variable("var2")
diff --git a/tests/test_solvers/test_scikits_solvers.py b/tests/test_solvers/test_scikits_solvers.py
index 6d0e37b8..59ab7d47 100644
--- a/tests/test_solvers/test_scikits_solvers.py
+++ b/tests/test_solvers/test_scikits_solvers.py
@@ -74,7 +74,7 @@ class TestScikitsSolver(unittest.TestCase):
         model.initial_conditions = {var: 1}
         model.initial_conditions_ydot = {var: 0.1}
         disc = model.default_discretisation
-        disc.process_model(model)
+        disc.process_model(model, model.default_geometry)
 
         # Solve
         solver = pybamm.ScikitsOdeSolver(tol=1e-9)
@@ -94,7 +94,7 @@ class TestScikitsSolver(unittest.TestCase):
         model.initial_conditions = {var1: 1, var2: 2}
         model.initial_conditions_ydot = {var1: 0.1, var2: 0.2}
         disc = model.default_discretisation
-        disc.process_model(model)
+        disc.process_model(model, model.default_geometry)
 
         # Solve
         solver = pybamm.ScikitsDaeSolver(tol=1e-8)
diff --git a/tests/test_solvers/test_scipy_solver.py b/tests/test_solvers/test_scipy_solver.py
index 599f3c84..b7d7960f 100644
--- a/tests/test_solvers/test_scipy_solver.py
+++ b/tests/test_solvers/test_scipy_solver.py
@@ -46,8 +46,10 @@ class TestScipySolver(unittest.TestCase):
 
         # create discretisation
         defaults = shared.TestDefaults1DMacro()
-        disc = shared.DiscretisationForTesting(defaults.mesh)
-        disc.process_model(model)
+        disc = shared.DiscretisationForTesting(
+            defaults.mesh_type, defaults.submesh_pts, defaults.submesh_types
+        )
+        disc.process_model(model, defaults.geometry)
         # Solve
         solver = pybamm.ScipySolver(tol=1e-8, method="RK45")
         t_eval = np.linspace(0, 1, 100)
diff --git a/tests/test_util.py b/tests/test_util.py
new file mode 100644
index 00000000..72d7034d
--- /dev/null
+++ b/tests/test_util.py
@@ -0,0 +1,62 @@
+#
+# Tests the Timer class.
+#
+# The code in this file is adapted from Pints
+# (see https://github.com/pints-team/pints)
+#
+import os
+import pybamm
+import unittest
+
+
+class TestUtil(unittest.TestCase):
+    """
+    Test the functionality in util.py
+    """
+
+    def test_load_function(self):
+        # Test filename ends in '.py'
+        with self.assertRaisesRegex(ValueError,
+                                    'Expected filename.py, but got doesnotendindotpy'):
+            pybamm.load_function('doesnotendindotpy')
+
+        # Test exception if absolute file not found
+        with self.assertRaisesRegex(ValueError,
+                                    'is an absolute path, but the file is not found'):
+            nonexistent_abs_file = os.path.join(os.getcwd(), 'i_dont_exist.py')
+            pybamm.load_function(nonexistent_abs_file)
+
+        # Test exception if relative file not found
+        with self.assertRaisesRegex(ValueError,
+                                    'cannot be found in the PyBaMM directory'):
+            pybamm.load_function('i_dont_exist.py')
+
+        # Test exception if relative file found more than once
+        with self.assertRaisesRegex(ValueError,
+                                    'found multiple times in the PyBaMM directory'):
+            pybamm.load_function('__init__.py')
+
+        # Test exception if no matching function found in module
+        with self.assertRaisesRegex(ValueError,
+                                    'No function .+ found in module .+'):
+            pybamm.load_function('process_symbol_bad_function.py')
+
+        # Test function load with absolute path
+        abs_test_path = os.path.join(os.getcwd(), 'tests', 'test_parameters', 'data',
+                                     'process_symbol_test_function.py')
+        self.assertTrue(os.path.isfile(abs_test_path))
+        func = pybamm.load_function(abs_test_path)
+        self.assertEqual(func(2), 246)
+
+        # Test function load with relative path
+        func = pybamm.load_function('process_symbol_test_function.py')
+        self.assertEqual(func(3), 369)
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    unittest.main()