void TestPottsMonolayerWithNonRandomSweep() throw (Exception)
    {
        EXIT_IF_PARALLEL;    // Potts simulations don't work in parallel because they depend on NodesOnlyMesh for writing.

        // Create a simple 2D PottsMesh
        PottsMeshGenerator<2> generator(6, 2, 2, 6, 2, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        MAKE_PTR(DifferentiatedCellProliferativeType, p_diff_type);
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasicRandom(cells, p_mesh->GetNumElements(), p_diff_type);

        // Create cell population
        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells);
        cell_population.SetUpdateNodesInRandomOrder(false);

        // Set up cell-based simulation
        OnLatticeSimulation<2> simulator(cell_population);
        simulator.SetOutputDirectory("TestSimplePottsMonolayerWithRandomSweep");
        simulator.SetEndTime(0.1);

        // Create update rules and pass to the simulation
        MAKE_PTR(VolumeConstraintPottsUpdateRule<2>, p_volume_constraint_update_rule);
        simulator.AddPottsUpdateRule(p_volume_constraint_update_rule);
        MAKE_PTR(AdhesionPottsUpdateRule<2>, p_adhesion_update_rule);
        simulator.AddPottsUpdateRule(p_adhesion_update_rule);

        // Run simulation
        TS_ASSERT_THROWS_NOTHING(simulator.Solve());
    }

    void TestGetLocationOfCellCentre() throw (Exception)
    {
        // Create a Potts-based cell population
        PottsMeshGenerator<2> generator(4, 2, 2, 4, 2, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells);

        // Test GetLocationOfCellCentre()
        AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();

        for (unsigned i=0; i<4; i++)
        {
            c_vector<double, 2> cell_location = cell_population.GetLocationOfCellCentre(*cell_iter);

            c_vector<double, 2> expected;
            expected(0) = 0.5 + 2*(i%2 != 0);
            expected(1) = 0.5 + 2*(i > 1);

            double drift = norm_2(cell_location-expected);
            TS_ASSERT_LESS_THAN(drift, 1e-6);

            ++cell_iter;
        }
    }

    void TestConstructor() throw(Exception)
    {
        EXIT_IF_PARALLEL;

        // Create a cell population
        HoneycombMeshGenerator generator(2, 2, 0);
        MutableMesh<2,2>* p_generating_mesh = generator.GetMesh();
        NodesOnlyMesh<2> mesh;
        mesh.ConstructNodesWithoutMesh(*p_generating_mesh, 1.5);

        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, mesh.GetNumNodes());

        NodeBasedCellPopulation<2> cell_population(mesh, cells);

        // Create a PDE handler object using this cell population
        CellBasedPdeHandler<2> pde_handler(&cell_population);

        // Test that member variables are initialised correctly
        TS_ASSERT_EQUALS(pde_handler.GetCellPopulation(), &cell_population);
        TS_ASSERT_EQUALS(pde_handler.GetWriteAverageRadialPdeSolution(), false);
        TS_ASSERT_EQUALS(pde_handler.GetWriteDailyAverageRadialPdeSolution(), false);
        TS_ASSERT_EQUALS(pde_handler.GetImposeBcsOnCoarseBoundary(), true);
        TS_ASSERT_EQUALS(pde_handler.GetNumRadialIntervals(), UNSIGNED_UNSET);
        TS_ASSERT(pde_handler.GetCoarsePdeMesh() == NULL);
    }

    /*
     * Note only solves one PDE
     */
    void TestMultipleCaBasedWithoutCoarseMeshUsingPdeHandlerOnCuboid() throw(Exception)
    {
        EXIT_IF_PARALLEL;

        // Create a simple 2D PottsMesh
        PottsMeshGenerator<2> generator(10, 0, 0, 10, 0, 0);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, 100);

        std::vector<unsigned> location_indices;
        for (unsigned i=0; i<100; i++)
        {
            location_indices.push_back(i);
        }

        // Create cell population
        MultipleCaBasedCellPopulation<2> cell_population(*p_mesh, cells, location_indices);

        // Set up cell-based simulation
        OnLatticeSimulation<2> simulator(cell_population);
        simulator.SetOutputDirectory("TestMultipleCaBasedCellPopulationWithPdesOnCuboid");
        simulator.SetDt(0.1);
        simulator.SetEndTime(1);

        // Set up a PDE with mixed boundary conditions (use zero uptake to check analytic solution)
        AveragedSourcePde<2> pde(cell_population, 0.0);
        ConstBoundaryCondition<2> bc(1.0);
        PdeAndBoundaryConditions<2> pde_and_bc(&pde, &bc, false);
        pde_and_bc.SetDependentVariableName("nutrient");

        // Pass this to the simulation object via a PDE handler
        CellBasedPdeHandlerOnCuboid<2> pde_handler(&cell_population);
        pde_handler.AddPdeAndBc(&pde_and_bc);
        pde_handler.SetImposeBcsOnCoarseBoundary(false);
        simulator.SetCellBasedPdeHandler(&pde_handler);

        // Solve the system
        simulator.Solve();

        // Test that PDE solver is working correctly
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            c_vector<double, 2> cell_location = simulator.rGetCellPopulation().GetLocationOfCellCentre(*cell_iter);

            if (cell_location[1] < 1e-6 || cell_location[1] > 9 - 1e-6)
            {
                TS_ASSERT_DELTA(cell_iter->GetCellData()->GetItem("nutrient"),1.0, 1e-2);
            }
            else
            {
                TS_ASSERT_LESS_THAN(1.0,cell_iter->GetCellData()->GetItem("nutrient"));
            }
        }
    }

    /*
     * The results may be visualized using {{{Visualize2dCentreCells}}} as described in the
     * previous test, with the results directory changed from {{{CellBasedDemo3}}} to {{{CellBasedDemo4}}}.
     *
     * == Test 5 - basic periodic mesh-based simulation ==
     *
     * We next show how to modify the previous test to implement a periodic boundary to the
     * left and right of the domain.
     */
    void TestMeshBasedMonolayerPeriodic() throw (Exception)
    {
        /* We now want to impose periodic boundaries on the domain. To do this we create a {{{Cylindrical2dMesh}}}
         * using a {{{CylindricalHoneycombMeshGenerator}}}.*/
        CylindricalHoneycombMeshGenerator generator(5, 2, 2); //**Changed**//
        Cylindrical2dMesh* p_mesh = generator.GetCylindricalMesh(); //**Changed**//

        /* Again we create one cell for each non ghost node. Note that we have changed back to using a {{{StochasticDurationCellCycleModel}}}.*/
        std::vector<unsigned> location_indices = generator.GetCellLocationIndices();
        std::vector<CellPtr> cells;
        MAKE_PTR(TransitCellProliferativeType, p_transit_type);
        CellsGenerator<StochasticDurationCellCycleModel, 2> cells_generator; //**Changed**//
        cells_generator.GenerateBasicRandom(cells, location_indices.size(), p_transit_type);

        /* We use the same {{{CellPopulation}}}, {{{CellBasedSimulation}}} (only changing the output directory and end time) and {{{Force}}} as before and run the simulation.*/
        MeshBasedCellPopulationWithGhostNodes<2> cell_population(*p_mesh, cells, location_indices);

        /* Again Paraview output is explicitly requested.*/
        cell_population.AddPopulationWriter<VoronoiDataWriter>();

        OffLatticeSimulation<2> simulator(cell_population);
        simulator.SetOutputDirectory("CellBasedDemo5"); //**Changed**//
        simulator.SetSamplingTimestepMultiple(12);
        simulator.SetEndTime(20.0); //**Changed**//

        MAKE_PTR(GeneralisedLinearSpringForce<2>, p_force);
        simulator.AddForce(p_force);

        simulator.Solve();

        /* The next two lines are for test purposes only and are not part of this tutorial.
         */
        TS_ASSERT_EQUALS(cell_population.GetNumRealCells(), 29u);
        TS_ASSERT_DELTA(SimulationTime::Instance()->GetTime(), 20.0, 1e-10);
    }

    void TestUpdateCellLocations()
    {
        // Create a simple 2D PottsMesh with two cells
        PottsMeshGenerator<2> generator(4, 2, 2, 2, 1, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        // Create cell population
        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells);

        // Set node selection to non-random lattice sweeping: this will loop over the nodes in index order
        TS_ASSERT_EQUALS(true, cell_population.GetUpdateNodesInRandomOrder());
        cell_population.SetUpdateNodesInRandomOrder(false);

        // Increase temperature: allows swaps to be more likely
        TS_ASSERT_EQUALS(cell_population.GetTemperature(),0.1);
        cell_population.SetTemperature(10.0);

        // Create a volume update rule and pass to the population
        MAKE_PTR(VolumeConstraintPottsUpdateRule<2>, p_volume_constraint_update_rule);
        cell_population.AddUpdateRule(p_volume_constraint_update_rule);

        // Commence lattice sweeping, updating where necessary
        // (Sensitive to changes in random number generation)
        cell_population.UpdateCellLocations(1.0);
        TS_ASSERT_EQUALS(cell_population.rGetCells().size(), 2u);
        TS_ASSERT_EQUALS(cell_population.rGetMesh().GetElement(0)->GetNumNodes(), 3u);
        TS_ASSERT_EQUALS(cell_population.rGetMesh().GetElement(1)->GetNumNodes(), 5u);
    }

    void TestUpdateCellLocationsRandomly()
    {
        // Create a simple 2D PottsMesh with two cells
        PottsMeshGenerator<2> generator(4, 2, 2, 2, 1, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        // Create cell population
        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells);

        // Set node selection to random lattice sweeping and (for coverage) random iteration over update rules
        cell_population.SetUpdateNodesInRandomOrder(true);
        cell_population.SetIterateRandomlyOverUpdateRuleCollection(true);

        // Increase temperature: allows swaps to be more likely
        TS_ASSERT_EQUALS(cell_population.GetTemperature(),0.1);
        cell_population.SetTemperature(10.0);

        // Create a volume update rule and pass to the population
        MAKE_PTR(VolumeConstraintPottsUpdateRule<2>, p_volume_constraint_update_rule);
        cell_population.AddUpdateRule(p_volume_constraint_update_rule);

        // Commence lattice sweeping, updating where necessary
        cell_population.UpdateCellLocations(1.0);

        // Note that these results differ to those in the above test due to extra random numbers being called
        // (Sensitive to changes in random number generation)
        TS_ASSERT_EQUALS(cell_population.rGetCells().size(), 2u);
        TS_ASSERT_EQUALS(cell_population.rGetMesh().GetElement(0)->GetNumNodes(), 4u);
        TS_ASSERT_EQUALS(cell_population.rGetMesh().GetElement(1)->GetNumNodes(), 4u);
    }

    void TestIsCellAssociatedWithADeletedLocation() throw (Exception)
    {
        // Create a Potts-based cell population but do not try to validate
        PottsMeshGenerator<2> generator(4, 2, 2, 4, 2, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();
        p_mesh->GetElement(0)->MarkAsDeleted();

        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells, false, false);

        // Test IsCellAssociatedWithADeletedLocation() method
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            bool is_deleted = cell_population.IsCellAssociatedWithADeletedLocation(*cell_iter);

            if (cell_population.GetLocationIndexUsingCell(*cell_iter) == 0)
            {
                TS_ASSERT_EQUALS(is_deleted, true);
            }
            else
            {
                TS_ASSERT_EQUALS(is_deleted, false);
            }
        }
    }

    void TestRemoveDeadCellsAndUpdate() throw(Exception)
    {
        // Create a simple 2D PottsMesh
        PottsMeshGenerator<2> generator(4, 2, 2, 4, 2, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        // Create cell population
        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells);

        // Test RemoveDeadCells() method
        TS_ASSERT_EQUALS(cell_population.GetNumElements(), 4u);

        cell_population.Begin()->Kill();

        TS_ASSERT_EQUALS(cell_population.RemoveDeadCells(), 1u);
        TS_ASSERT_EQUALS(cell_population.GetNumElements(), 3u);
        TS_ASSERT_EQUALS(p_mesh->GetNumElements(), 3u);
        TS_ASSERT_EQUALS(p_mesh->GetNumAllElements(), 4u);

        // Test that Update() throws no errors
        TS_ASSERT_THROWS_NOTHING(cell_population.Update());
    }

    void TestChemotaxisPottsUpdateRuleIn2d() throw (Exception)
    {
        // Create a simple 2D PottsMesh with 2 elements
        PottsMeshGenerator<2> generator(4, 1, 2, 4, 1, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        // Create cell population
        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells);

        // Create an update law system
        ChemotaxisPottsUpdateRule<2> chemotaxis_update;

        // Test EvaluateHamiltonianContribution()

        // target site above current site
        double contribution = chemotaxis_update.EvaluateHamiltonianContribution(10, 14, cell_population);
        TS_ASSERT_DELTA(contribution, -0.2, 1e-6);

        // target site below current site
        contribution = chemotaxis_update.EvaluateHamiltonianContribution(6, 2, cell_population);
        TS_ASSERT_DELTA(contribution, 0.2, 1e-6);

        // target site diagonally above current site (to right)
        contribution = chemotaxis_update.EvaluateHamiltonianContribution(10, 15, cell_population);
        TS_ASSERT_DELTA(contribution, -0.4, 1e-6);

        // target site diagonally above current site (to left)
        contribution = chemotaxis_update.EvaluateHamiltonianContribution(10, 13, cell_population);
        TS_ASSERT_DELTA(contribution, 0.0, 1e-6);
    }

    void TestCellwiseDataGradientVerySmallMesh() throw(Exception)
    {
        // Create a simple mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_2_elements");
        MutableMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        // Create a cell population
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, mesh.GetNumNodes());
        MeshBasedCellPopulation<2> cell_population(mesh, cells);

        // Set up data: C(x,y) = x^2
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            double x = mesh.GetNode(i)->rGetLocation()[0];
            CellPtr p_cell = cell_population.GetCellUsingLocationIndex(mesh.GetNode(i)->GetIndex());
            p_cell->GetCellData()->SetItem("x^2", x*x);
        }

        CellwiseDataGradient<2> gradient;
        gradient.SetupGradients(cell_population, "x^2");

        // With the algorithm being used, the numerical gradient is (1,0)
        // for each of the nodes
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            TS_ASSERT_DELTA(gradient.rGetGradient(i)(0), 1.0, 1e-9);
            TS_ASSERT_DELTA(gradient.rGetGradient(i)(1), 0.0, 1e-9);
        }
    }

    /*
     * EMPTYLINE
     *
     * To visualize the results, open a new terminal, {{{cd}}} to the Chaste directory,
     * then {{{cd}}} to {{{anim}}}. Then do: {{{java Visualize2dCentreCells /tmp/$USER/testoutput/NodeBasedMonolayer/results_from_time_0}}}.
     * we need to select the 'Cells as circles` option to be able to visualize the cells, as opposed
     * to just the centres.
     * We may have to do: {{{javac Visualize2dCentreCells.java}}} beforehand to create the
     * java executable.
     *
     * Alternatively, to view in Paraview, load the file {{{/tmp/$USER/testoutput/NodeBasedMonolayer/results_from_time_0/results.pvd}}}
     * and add glyphs to represent cells. An option is to use 3D spherical glyphs and then make a planar cut.
     * Note that, for larger simulations, you may need to unclick "Mask Points" (or similar) so as not to limit the number of glyphs
     * displayed by Paraview.
     *
     *
     *
     * EMPTYLINE
     *
     * == Test 2 - a basic node-based simulation in 3D ==
     *
     * EMPTYLINE
     *
     * In the second test we run a simple node-based simulation in 3D. This is very similar
     * to the 2D test with the dimension template (<2,2> and <2>) changed from 2 to 3 and instead of using a mesh
     * generator we generate the nodes directly.
     */
    void TestSpheroid()
    {
        /** The next line is needed because we cannot currently run node based simulations in parallel. */
        EXIT_IF_PARALLEL;

        /*
         * First, we generate a nodes only mesh. This time we specify the nodes manually by first
         * creating a vector of nodes. */
        std::vector<Node<3>*> nodes;
        /* We then create some nodes to add to this vector. */
        nodes.push_back(new Node<3>(0u,  false,  0.5, 0.0, 0.0));
        nodes.push_back(new Node<3>(1u,  false,  -0.5, 0.0, 0.0));
        nodes.push_back(new Node<3>(2u,  false,  0.0, 0.5, 0.0));
        nodes.push_back(new Node<3>(3u,  false,  0.0, -0.5, 0.0));
        /* Finally a {{{NodesOnlyMesh}}} is created and the vector of nodes is passed to
         * the {{{ConstructNodesWithoutMesh}}} method. */
        NodesOnlyMesh<3> mesh;
        /* To run node-based simulations you need to define a cut off length (second argument in
         * {{{ConstructNodesWithoutMesh}}}), which defines the connectivity of the nodes by defining
         * a radius of interaction. */
        mesh.ConstructNodesWithoutMesh(nodes, 1.5);

        /*
         * Having created a mesh, we now create a {{{std::vector}}} of {{{CellPtr}}}s.
         * As before, we do this with the `CellsGenerator` helper class (this time with dimension 3).
         */
        std::vector<CellPtr> cells;
        MAKE_PTR(TransitCellProliferativeType, p_transit_type);
        CellsGenerator<UniformCellCycleModel, 3> cells_generator;
        cells_generator.GenerateBasicRandom(cells, mesh.GetNumNodes(), p_transit_type);

        /* We make a {{{NodeBasedCellPopulation}}} (this time with dimension 3) as before.
         */
        NodeBasedCellPopulation<3> cell_population(mesh, cells);

        /* We then pass in the cell population into an {{{OffLatticeSimulation}}},
         * (this time with dimension 3) and set the output directory, output multiple and end time. */
        OffLatticeSimulation<3> simulator(cell_population);
        simulator.SetOutputDirectory("NodeBasedSpheroid");
        simulator.SetSamplingTimestepMultiple(12);
        simulator.SetEndTime(10.0);

        /* Again we create a force law (this time with dimension 3), and pass it to the {{{OffLatticeSimulation}}}.*/
        MAKE_PTR(GeneralisedLinearSpringForce<3>, p_force);
        simulator.AddForce(p_force);

        /* To run the simulation, we call {{{Solve()}}}. */
        simulator.Solve();

        /* The next two lines are for test purposes only and are not part of this tutorial.
         */
        TS_ASSERT_EQUALS(cell_population.GetNumRealCells(), 8u);
        TS_ASSERT_DELTA(SimulationTime::Instance()->GetTime(), 10.0, 1e-10);

        /* To avoid memory leaks, we conclude by deleting any pointers that we created in the test.*/
        for (unsigned i=0; i<nodes.size(); i++)
        {
            delete nodes[i];
        }
    }

    /*
     * == Using the SRN model in a cell-based simulation ==
     *
     * We conclude with a brief test demonstrating how {{{MySrnModel}}} can be used
     * in a cell-based simulation.
     */
    void TestOffLatticeSimulationWithMySrnModel() throw(Exception)
    {
        /* We use the honeycomb vertex mesh generator to create a vertex mesh.
         */
        HoneycombVertexMeshGenerator generator(2, 2);
        MutableVertexMesh<2,2>* p_mesh = generator.GetMesh();

        /* Next, we create some cells. First, define the cells vector. */
        std::vector<CellPtr> cells;
        /* We must create a shared_ptr to a {{{CellMutationState}}} with which to bestow the cells.
         * We make use of the macro MAKE_PTR to do this: the first argument is the class and
         * the second argument is the name of the shared_ptr. */
        MAKE_PTR(WildTypeCellMutationState, p_state);
        MAKE_PTR(StemCellProliferativeType, p_stem_type);
        /* Then we loop over the nodes. */
        for (unsigned i=0; i<p_mesh->GetNumElements(); i++)
        {
            /* For each node we create a cell with our SRN model and simple Stochastic cell cycle model. */
            StochasticDurationCellCycleModel* p_cell_cycle_model = new StochasticDurationCellCycleModel();
            MySrnModel* p_srn_model = new MySrnModel;

            /* We choose to initialise the concentrations to random levels in each cell. */
            std::vector<double> initial_conditions;
            initial_conditions.push_back(1.0-2.0*RandomNumberGenerator::Instance()->ranf());
            initial_conditions.push_back(1.0-2.0*RandomNumberGenerator::Instance()->ranf());
            p_srn_model->SetInitialConditions(initial_conditions);

            CellPtr p_cell(new Cell(p_state, p_cell_cycle_model, p_srn_model));
            p_cell->SetCellProliferativeType(p_stem_type);


            /* Now, we define a random birth time, chosen from [-T,0], where
             * T = t,,1,, + t,,2,,, where t,,1,, is a parameter representing the G,,1,, duration
             * of a stem cell, and t,,2,, is the basic S+G,,2,,+M phases duration.
             */
            double birth_time = - RandomNumberGenerator::Instance()->ranf() * (p_cell_cycle_model->GetStemCellG1Duration() + p_cell_cycle_model->GetSG2MDuration());
            /* We then set the birth time and push the cell back into the vector of cells. */
            p_cell->SetBirthTime(birth_time);
            cells.push_back(p_cell);
        }


        /* Now that we have defined the mesh and cells, we can define the cell population, forces, areas modifier, and simulation
         * in the same way as the other tutorials. */
        VertexBasedCellPopulation<2> cell_population(*p_mesh, cells);

        OffLatticeSimulation<2> simulator(cell_population);
        simulator.SetOutputDirectory("TestOffLatticeSimulationWithMySrnModel");
        simulator.SetEndTime(10.0);
        simulator.SetSamplingTimestepMultiple(50);

        MAKE_PTR(NagaiHondaForce<2>, p_force);
        simulator.AddForce(p_force);

        MAKE_PTR(SimpleTargetAreaModifier<2>, p_growth_modifier);
        simulator.AddSimulationModifier(p_growth_modifier);

        /* Finally to run the simulation, we call {{{Solve()}}}. */
        simulator.Solve();
    }

    void TestAllCases()
    {
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_2_elements");
        MutableMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        TS_ASSERT_DELTA(mesh.GetAngleBetweenNodes(2,0), -0.75*M_PI, 1e-12);
        TS_ASSERT_DELTA(mesh.GetAngleBetweenNodes(2,1), -0.5*M_PI, 1e-12);

        CellsGenerator<FixedG1GenerationalCellCycleModel, 2> cells_generator;
        std::vector<CellPtr> cells;
        cells_generator.GenerateBasic(cells, mesh.GetNumNodes());

        MeshBasedCellPopulation<2> cell_population(mesh, cells);

        MAKE_PTR(GeneralisedLinearSpringForce<2>, p_force);
        std::vector<boost::shared_ptr<AbstractTwoBodyInteractionForce<2> > > force_collection;
        force_collection.push_back(p_force);

        DiscreteSystemForceCalculator calculator(cell_population, force_collection);

        double epsilon = 0.5*M_PI;
        calculator.SetEpsilon(epsilon);

        TS_ASSERT_THROWS_NOTHING(calculator.GetSamplingAngles(2));
    }

    void TestAddingUpdateRules() throw(Exception)
    {
        // Create a simple 2D PottsMesh with one cell
        PottsMeshGenerator<2> generator(2, 1, 2, 2, 1, 2);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        // Create cell population
        PottsBasedCellPopulation<2> cell_population(*p_mesh, cells);

        // Test we have the correct number of cells and elements
        TS_ASSERT_EQUALS(cell_population.GetNumElements(), 1u);
        TS_ASSERT_EQUALS(cell_population.rGetCells().size(), 1u);

        // Create an update rule and pass to the population
        MAKE_PTR(VolumeConstraintPottsUpdateRule<2>, p_volume_constraint_update_rule);
        cell_population.AddUpdateRule(p_volume_constraint_update_rule);

        // Check the update rules are correct
        std::vector<boost::shared_ptr<AbstractPottsUpdateRule<2> > > update_rule_collection = cell_population.rGetUpdateRuleCollection();
        TS_ASSERT_EQUALS(update_rule_collection.size(),1u);
        TS_ASSERT_EQUALS((*update_rule_collection[0]).GetIdentifier(), "VolumeConstraintPottsUpdateRule-2");
    }

    void TestCellPopulationIteratorWithNoCells() throw(Exception)
    {
        EXIT_IF_PARALLEL;    // This test doesn't work in parallel.

        SimulationTime* p_simulation_time = SimulationTime::Instance();
        p_simulation_time->SetEndTimeAndNumberOfTimeSteps(10.0, 1);

        // Create a simple mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_128_elements");
        TetrahedralMesh<2,2> generating_mesh;
        generating_mesh.ConstructFromMeshReader(mesh_reader);

        // Convert this to a NodesOnlyMesh
        NodesOnlyMesh<2> mesh;
        mesh.ConstructNodesWithoutMesh(generating_mesh, 1.2);

        // Create vector of cell location indices
        std::vector<unsigned> cell_location_indices;
        cell_location_indices.push_back(80);

        // Create a single cell
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, cell_location_indices.size());
        cells[0]->StartApoptosis();

        // Create a cell population
        NodeBasedCellPopulationWithParticles<2> cell_population(mesh, cells, cell_location_indices);

        // Iterate over cell population and check there is a single cell
        unsigned counter = 0;
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
                cell_iter != cell_population.End();
                ++cell_iter)
        {
            counter++;
        }
        TS_ASSERT_EQUALS(counter, 1u);
        TS_ASSERT_EQUALS(cell_population.rGetCells().empty(), false);

        // Increment simulation time and update cell population
        p_simulation_time->IncrementTimeOneStep();

        unsigned num_cells_removed = cell_population.RemoveDeadCells();
        TS_ASSERT_EQUALS(num_cells_removed, 1u);

        cell_population.Update();

        // Iterate over cell population and check there are now no cells
        counter = 0;
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
                cell_iter != cell_population.End();
                ++cell_iter)
        {
            counter++;
        }
        TS_ASSERT_EQUALS(counter, 0u);
        TS_ASSERT_EQUALS(cell_population.rGetCells().empty(), true);
    }

    /*
     * Cellular birth has been tested in TestCaSingleCellWithBirth for one cell per lattice site.
     * This test adds to the above by further testing cellular birth considering multiple cells per lattice site.
     * A  two-lattice mesh was created and only one lattice had free space to add one daughter cell.
     */
    void TestMultipleCellsPerLatticeSiteWithBirth() throw (Exception)
    {
        EXIT_IF_PARALLEL;

        // Create a simple 2D PottsMesh
        PottsMeshGenerator<2> generator(2, 0, 0, 1, 0, 0);
        PottsMesh<2>* p_mesh = generator.GetMesh();

        // Create cells
        std::vector<CellPtr> cells;
        MAKE_PTR(StemCellProliferativeType, p_stem_type);
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasicRandom(cells, 3, p_stem_type);

        // Specify where cells lie
        std::vector<unsigned> location_indices;
        location_indices.push_back(0);
        location_indices.push_back(0);
        location_indices.push_back(1);

        // Create cell population
        CaBasedCellPopulation<2> cell_population(*p_mesh, cells, location_indices, 2);

        // Set up cell-based simulation
        OnLatticeSimulation<2> simulator(cell_population);
        std::string output_directory = "TestMultipleCellsPerLatticeSiteWithBirth";
        simulator.SetOutputDirectory(output_directory);
        simulator.SetDt(0.1);
        simulator.SetEndTime(40);

        // Add update rule
        MAKE_PTR(DiffusionCaUpdateRule<2u>, p_diffusion_update_rule);
        p_diffusion_update_rule->SetDiffusionParameter(0.5);
        simulator.AddCaUpdateRule(p_diffusion_update_rule);

        // Run simulation
        simulator.Solve();

        // Check the number of cells
        TS_ASSERT_EQUALS(simulator.rGetCellPopulation().GetNumRealCells(), 4u);

        // Test no deaths and some births
        TS_ASSERT_EQUALS(simulator.GetNumBirths(), 1u);
        TS_ASSERT_EQUALS(simulator.GetNumDeaths(), 0u);

#ifdef CHASTE_VTK
        // Test that VTK writer has produced a file
        OutputFileHandler output_file_handler(output_directory, false);
        std::string results_dir = output_file_handler.GetOutputDirectoryFullPath();

        // Initial condition file
        FileFinder vtk_file(results_dir + "results_from_time_0/results_0.vtu", RelativeTo::Absolute);
        TS_ASSERT(vtk_file.Exists());

        // Final file
        FileFinder vtk_file2(results_dir + "results_from_time_0/results_400.vtu", RelativeTo::Absolute);
        TS_ASSERT(vtk_file2.Exists());
#endif //CHASTE_VTK
    }

    void TestAdhesionPottsUpdateRuleIn3d() throw (Exception)
    {
        // Create a simple 2D PottsMesh with 2 elements
        PottsMeshGenerator<3> generator(4, 2, 2, 2, 1, 2, 4, 1, 2, true);
        PottsMesh<3>* p_mesh = generator.GetMesh();

        TS_ASSERT_EQUALS(p_mesh->GetNumElements(), 2u);
        TS_ASSERT_EQUALS(p_mesh->GetNumNodes(), 32u);

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 3> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        // Label cells 0 and 1
        MAKE_PTR(CellLabel, p_label);
        cells[0]->AddCellProperty(p_label);
        cells[1]->AddCellProperty(p_label);

        // Create cell population
        PottsBasedCellPopulation<3> cell_population(*p_mesh, cells);

        // Create an update law system
        AdhesionPottsUpdateRule<3> adhesion_update;

        // Set the parameters to facilitate calculations below.
        adhesion_update.SetCellCellAdhesionEnergyParameter(0.1);
        adhesion_update.SetCellBoundaryAdhesionEnergyParameter(0.2);

        // Test EvaluateHamiltonianContribution(). Note that the boundary of the domain is a void not medium.

        double gamma_cell_cell = adhesion_update.GetCellCellAdhesionEnergyParameter();
        double gamma_cell_boundary = adhesion_update.GetCellBoundaryAdhesionEnergyParameter();

        // Both points lie within cell 0
        TS_ASSERT_THROWS_THIS(adhesion_update.EvaluateHamiltonianContribution(0, 1, cell_population),
                              "The current node and target node must not be in the same element.");

        // Both points lie within cell 1
        TS_ASSERT_THROWS_THIS(adhesion_update.EvaluateHamiltonianContribution(2, 3, cell_population),
                              "The current node and target node must not be in the same element.");

        // Both points lie within cell medium
        TS_ASSERT_THROWS_THIS(adhesion_update.EvaluateHamiltonianContribution(16, 17, cell_population),
                               "At least one of the current node or target node must be in an element.");

        // Current site in cell 0; target site in cell medium
        double contribution = adhesion_update.EvaluateHamiltonianContribution(9, 17, cell_population);
        TS_ASSERT_DELTA(contribution, 3.0*gamma_cell_boundary, 1e-6);

        // Current site in cell medium; target site in cell 0
        contribution = adhesion_update.EvaluateHamiltonianContribution(17, 9, cell_population);
        TS_ASSERT_DELTA(contribution, 3.0*gamma_cell_boundary - gamma_cell_cell, 1e-6);

        // Current site in cell 0; target site in cell 1
        contribution = adhesion_update.EvaluateHamiltonianContribution(9, 10, cell_population);
        TS_ASSERT_DELTA(contribution, 2.0*gamma_cell_cell, 1e-6);
    }

    void TestSurfaceAreaConstraintPottsUpdateRuleIn3d() throw (Exception)
    {
        // Create a simple 3D PottsMesh with 2 elements
        PottsMeshGenerator<3> generator(4, 2, 2, 4, 1, 2, 4, 1, 2, true);
        PottsMesh<3>* p_mesh = generator.GetMesh();

        TS_ASSERT_EQUALS(p_mesh->GetNumElements(), 2u);
        TS_ASSERT_EQUALS(p_mesh->GetNumNodes(), 64u);

        // Create cells
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 3> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumElements());

        // Create cell population
        PottsBasedCellPopulation<3> cell_population(*p_mesh, cells);

        // Create an update law system
        SurfaceAreaConstraintPottsUpdateRule<3> surface_area_constraint;

        surface_area_constraint.SetDeformationEnergyParameter(0.5);
        surface_area_constraint.SetMatureCellTargetSurfaceArea(24.0);

        // Test EvaluateHamiltonianContribution()
        double gamma = surface_area_constraint.GetDeformationEnergyParameter();

        // Both points lie within cell 0
        TS_ASSERT_THROWS_THIS(surface_area_constraint.EvaluateHamiltonianContribution(0, 1, cell_population),
                              "The current node and target node must not be in the same element.");

        // Both points lie within cell 1
        TS_ASSERT_THROWS_THIS(surface_area_constraint.EvaluateHamiltonianContribution(2, 3, cell_population),
                              "The current node and