As with any contribution to PyBaMM, please follow the workflow in CONTRIBUTING.md. In particular, start by creating an issue to discuss what you want to do - this is a good way to avoid wasted coding hours!

All models in PyBaMM are implemented as expression trees.
After it has been created and parameters have been set, the model is passed to the `pybamm.Discretisation`

class,
which converts it into a linear algebra form.
For example, the object:

```
grad(u)
```

might get converted to a Matrix-Vector multiplication:

```
Matrix(100,100) @ y[0:100]
```

(in Python 3.5+, @ means matrix multiplication, while * is elementwise product).
The `pybamm.Discretisation`

class is a wrapper that iterates through the different parts of the model,
performing the trivial conversions (e.g. Addition –> Addition),
and calls upon spatial methods to perform the harder conversions (e.g. grad(u) –> Matrix * StateVector, SpatialVariable –> Vector, etc).

Hence SpatialMethod classes only need to worry about the specific conversions, and `pybamm.Discretisation`

deals with the rest.

To add a new spatial method (e.g. My Fast Method), first create a new file (`my_fast_method.py`

) in `pybamm/spatial_methods/`

,
with a single class that inherits from `pybamm.SpatialMethod`

, such as:

```
class MyFastMethod(pybamm.SpatialMethod):
```

and add the class to `pybamm/__init__.py`

:

```
from .spatial_methods.my_fast_method import MyFastMethod
```

You can then start implementing the spatial method by adding functions to the class.
In particular, any spatial method *must* have the following functions (from the base class `pybamm.SpatialMethod`

):

`pybamm.SpatialMethod.gradient()`

`pybamm.SpatialMethod.divergence()`

`pybamm.SpatialMethod.integral()`

`pybamm.SpatialMethod.indefinite integral()`

`pybamm.SpatialMethod.boundary_value_or_flux()`

Optionally, a new spatial method can also overwrite the default behaviour for the following functions:

`pybamm.SpatialMethod.spatial_variable()`

`pybamm.SpatialMethod.broadcast()`

`pybamm.SpatialMethod.mass_matrix()`

`pybamm.SpatialMethod.process_binary_operators()`

`pybamm.SpatialMethod.concatenation()`

For an example of an existing spatial method implementation, see the Finite Volume API docs and notebook.

For the new spatial method to be added to PyBaMM, you must add unit tests to demonstrate that it behaves as expected
(see, for example, the Finite Volume unit tests).
The best way to get started would be to create a file `test_my_fast_method.py`

in `tests/unit/test_spatial_methods/`

that performs at least the
following checks:

- Operations return objects that have the expected shape
- Standard operations behave as expected, e.g. (in 1D) grad(x^2) = 2*x, integral(sin(x), 0, pi) = 2
- (more advanced) make sure that the operations converge at the correct rate to known analytical solutions as you decrease the grid size

In theory, any existing model can now be discretised using `MyFastMethod`

instead of their default spatial methods, with no extra work from here.
To test this, add something like the following test to one of the model test files
(e.g. DFN):

```
def test_my_fast_method(self):
model = pybamm.lithium_ion.DFN()
spatial_methods = {
"macroscale": pybamm.MyFastMethod,
"negative particle": pybamm.MyFastMethod,
"positive particle": pybamm.MyFastMethod,
}
modeltest = tests.StandardModelTest(model, spatial_methods=spatial_methods)
modeltest.test_all()
```

This will check that the model can run with the new spatial method (but not that it gives a sensible answer!).

Once you have performed the above checks, you are almost ready to merge your code into the core PyBaMM - see CONTRIBUTING.md workflow for how to do this.