Adding a custom model
=====================
If :ref:`adding a custom boundary condition` is insufficient for your needs, it is also relatively straightforward to implement a custom model.
In this case, however, the :class:`.Discretization` class will need to be modified (or the appropriate methods needs to be replaced in an instance of the class).
The methods that most likely need to be changed are :meth:`linear_part`, :meth:`nonlinear_part`, :meth:`mass_matrix` and :meth:`boundaries`, or functions called from any of these functions.

As an example, the 2D nondimensional incompressible Stokes equations of the form

.. math::
   \frac{1}{\mathrm{Re}} \nabla^2 u - \nabla p &= 0\\
   \nabla \cdot u &= 0

could be implemented as

.. code-block:: Python

    def linear_part(self):
        Re = self.get_parameter('Reynolds Number', 1.0)
        return 1 / Re * (self.u_xx() + self.u_yy()
                         + self.v_xx() + self.v_yy()) \
            - (self.p_x() + self.p_y()) \
            + self.div()

    def nonlinear_part(self, state):
        atom = numpy.zeros((self.nx, self.ny, self.nz, self.dof, self.dof, 3, 3, 3))
        return atom, atom

We could then add a moving lid to the top using

.. code-block:: Python

    def boundaries(self, atom):
        boundary_conditions = BoundaryConditions(
            self.nx, self.ny, self.nz, self.dim, self.dof, self.x, self.y, self.z)

        lid_velocity = self.parameters['Lid Velocity']
        boundary_conditions.moving_lid_north(atom, lid_velocity)

        boundary_conditions.no_slip_east(atom)
        boundary_conditions.no_slip_west(atom)

        boundary_conditions.no_slip_south(atom)

        return boundary_conditions.get_forcing()

..
    Explicitly enable math mode
.. math::