C++ MutableMesh::ConstructFromMeshReader方法代码示例

    void TestDeleteFirstOrder() throw(Exception)
    {
#ifdef CHASTE_VTK
        VtkMeshReader<1,3> mesh_reader("lung/test/data/TestSubject002MajorAirways.vtu");

        MutableMesh<1,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        //Assign valid radii
        for(unsigned node_index = 0; node_index < mesh.GetNumNodes(); ++node_index)
        {
            mesh.GetNode(node_index)->AddNodeAttribute(1.0);
        }

        TS_ASSERT_EQUALS(mesh.GetNumNodes(), 12065u);
        TS_ASSERT_EQUALS(mesh.GetNumElements(), 12064u);

        MajorAirwaysCentreLinesCleaner cleaner(mesh, 0u);
        cleaner.CleanUsingHorsfieldOrder(1u);

        NodeMap node_map(mesh.GetNumAllNodes());
        mesh.ReIndex(node_map);

        //Confirmed visually to be correct
        TS_ASSERT_EQUALS(mesh.GetNumNodes(), 3683u);
        TS_ASSERT_EQUALS(mesh.GetNumElements(), 3682u);

// Uncomment to visualise
//        VtkMeshWriter<1,3> mesh_writer("TestMajorAirwaysCentreLinesCleaner", "Novartis002Trimmed");
//        mesh_writer.WriteFilesUsingMesh(mesh);
#endif //CHASTE_VTK
    }

    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);
        }
    }

    void TestDeleteOrderSimpleMesh() throw(Exception)
    {
        TrianglesMeshReader<1,3> mesh_reader("lung/test/airway_generation/data/test_major_airways_mesh");
        MutableMesh<1,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        //Assign valid radii
        for(unsigned node_index = 0; node_index < mesh.GetNumNodes(); ++node_index)
        {
            mesh.GetNode(node_index)->rGetNodeAttributes()[0] = 1.0;
        }

        TS_ASSERT_EQUALS(mesh.GetNumNodes(), 6u);
        TS_ASSERT_EQUALS(mesh.GetNumElements(), 5u);

        MajorAirwaysCentreLinesCleaner cleaner(mesh, 0u);
        cleaner.CleanUsingHorsfieldOrder(1u);

        NodeMap node_map(mesh.GetNumAllNodes());
        mesh.ReIndex(node_map);

        TS_ASSERT_EQUALS(mesh.GetNumNodes(), 2u);
        TS_ASSERT_EQUALS(mesh.GetNumElements(), 1u);
        TS_ASSERT_DELTA(mesh.GetNode(0u)->rGetLocation()[0], 0.0, 1e-6);
        TS_ASSERT_DELTA(mesh.GetNode(0u)->rGetLocation()[1], 0.0, 1e-6);
        TS_ASSERT_DELTA(mesh.GetNode(0u)->rGetLocation()[2], -2.0, 1e-6);
        TS_ASSERT_DELTA(mesh.GetNode(1u)->rGetLocation()[0], 0.0, 1e-6);
        TS_ASSERT_DELTA(mesh.GetNode(1u)->rGetLocation()[1], 0.0, 1e-6);
        TS_ASSERT_DELTA(mesh.GetNode(1u)->rGetLocation()[2], 0.0, 1e-6);
    }

    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 TestSloughingCellKillerTopAndSides() throw(Exception)
    {
        // Create mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_128_elements");
        MutableMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);
        mesh.Translate(-0.25,-0.25);

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

        // Create cell population
        MeshBasedCellPopulation<2> cell_population(mesh, cells);

        // Create cell killer and kill cells
        SloughingCellKiller<2> sloughing_cell_killer(&cell_population, 0.5, true, 0.5);

        sloughing_cell_killer.CheckAndLabelCellsForApoptosisOrDeath();

        TS_ASSERT_EQUALS(sloughing_cell_killer.GetIdentifier(), "SloughingCellKiller-2");

        // Check that cells were labelled for death correctly
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            double x = cell_population.GetLocationOfCellCentre(*cell_iter)[0];
            double y = cell_population.GetLocationOfCellCentre(*cell_iter)[1];

            if ((x<0) || (x>0.5) || (y>0.5))
            {
                TS_ASSERT_EQUALS(cell_iter->IsDead(), true);
            }
            else
            {
                TS_ASSERT_EQUALS(cell_iter->IsDead(), false);
            }
        }

        cell_population.RemoveDeadCells();

        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            double x = cell_population.GetLocationOfCellCentre(*cell_iter)[0];
            double y = cell_population.GetLocationOfCellCentre(*cell_iter)[1];

            TS_ASSERT_LESS_THAN_EQUALS(x, 0.5);
            TS_ASSERT_LESS_THAN_EQUALS(y, 0.5);
        }
    }

    void TestUseInPopulationWriteResultsToFile()
    {
        EXIT_IF_PARALLEL;

        // Set up SimulationTime (needed if VTK is used)
        SimulationTime::Instance()->SetEndTimeAndNumberOfTimeSteps(1.0, 1);

        // Create a simple mesh-based cell population, comprising various cell types in various cell cycle phases
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_4_elements");
        MutableMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

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

        MeshBasedCellPopulation<2> cell_population(mesh, cells);

        // Create an instance of a Wnt concentration
        WntConcentration<2>::Instance()->SetType(LINEAR);
        WntConcentration<2>::Instance()->SetCellPopulation(cell_population);
        WntConcentration<2>::Instance()->SetCryptLength(1.0);

        cell_population.InitialiseCells();

        // This is where we add the writer
        cell_population.AddCellWriter<CellBetaCateninWriter>();

        // This method is usually called by Update()
        cell_population.CreateVoronoiTessellation();

        std::string output_directory = "TestUseInPopulationWriteResultsToFile";
        OutputFileHandler output_file_handler(output_directory, false);

        cell_population.OpenWritersFiles(output_file_handler);
        cell_population.WriteResultsToFiles(output_directory);

        SimulationTime::Instance()->IncrementTimeOneStep();
        cell_population.Update();

        cell_population.WriteResultsToFiles(output_directory);
        cell_population.CloseWritersFiles();

        // Compare output with saved file of what they should look like
        std::string results_dir = output_file_handler.GetOutputDirectoryFullPath();

        FileComparison comparer(results_dir + "results.vizbetacatenin","crypt/test/data/TestCellBetaCateninWriter/results2.vizbetacatenin");
        TS_ASSERT(comparer.CompareFiles());
    }

    void TestRemoveIsolatedNodesSimpleMesh() throw(Exception)
    {
        TrianglesMeshReader<1,3> mesh_reader("lung/test/airway_generation/data/test_isolated_nodes_major_airways_mesh");
        MutableMesh<1,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        //Assign valid radii
        for(unsigned node_index = 0; node_index < mesh.GetNumNodes(); ++node_index)
        {
          mesh.GetNode(node_index)->rGetNodeAttributes()[0] = 1.0;
        }

        TS_ASSERT_EQUALS(mesh.GetNumNodes(), 6u);
        TS_ASSERT_EQUALS(mesh.GetNumElements(), 3u);

        MajorAirwaysCentreLinesCleaner cleaner(mesh, 0u);
        cleaner.CleanIsolatedNodes();

        TS_ASSERT_EQUALS(mesh.GetNumNodes(), 4u);
        TS_ASSERT_EQUALS(mesh.GetNumElements(), 3u);
    }

    void TestSolvePdeAndWriteResultsToFileCoarsePdeMeshNeumann() throw(Exception)
    {
        EXIT_IF_PARALLEL;

        // Set up SimulationTime
        SimulationTime::Instance()->SetEndTimeAndNumberOfTimeSteps(0.05, 6);

        // Create a cigar-shaped mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/disk_522_elements");
        MutableMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);
        mesh.Scale(5.0, 1.0);

        // 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);

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

        // Set up PDE and pass to handler
        AveragedSourcePde<2> pde(cell_population, -0.01);
        ConstBoundaryCondition<2> bc(0.0);
        PdeAndBoundaryConditions<2> pde_and_bc(&pde, &bc, true); // Last boolean specifies Neuman conditions
        pde_and_bc.SetDependentVariableName("variable");

        pde_handler.AddPdeAndBc(&pde_and_bc);

        // Solve PDEs on a coarse mesh
        ChastePoint<2> lower(0.0, 0.0);
        ChastePoint<2> upper(50.0, 50.0);
        ChasteCuboid<2> cuboid(lower, upper);
        pde_handler.UseCoarsePdeMesh(10.0, cuboid, true);
        pde_handler.SetImposeBcsOnCoarseBoundary(false);

        // For coverage, provide an initial guess for the solution
        std::vector<double> data(pde_handler.mpCoarsePdeMesh->GetNumNodes());
        for (unsigned i=0; i<pde_handler.mpCoarsePdeMesh->GetNumNodes(); i++)
        {
            data[i] = 1.0;
        }

        Vec vector = PetscTools::CreateVec(data);
        pde_and_bc.SetSolution(vector);

        // Open result file ourselves
        OutputFileHandler output_file_handler("TestWritePdeSolution", false);
        pde_handler.mpVizPdeSolutionResultsFile = output_file_handler.OpenOutputFile("results.vizpdesolution");

        // Solve PDE (set sampling timestep multiple to be large doesn't do anything as always output on 1st timestep)
        pde_handler.SolvePdeAndWriteResultsToFile(10);

        // Close result file ourselves
        pde_handler.mpVizPdeSolutionResultsFile->close();

        // Test that boundary cells experience the right boundary condition
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            if (cell_population.GetNodeCorrespondingToCell(*cell_iter)->IsBoundaryNode())
            {
                TS_ASSERT_DELTA(cell_iter->GetCellData()->GetItem("variable"), 0.0, 1e-1);
            }
        }
    }

    void TestSolvePdeAndWriteResultsToFileWithoutCoarsePdeMeshNeumann() throw(Exception)
    {
        EXIT_IF_PARALLEL;

        // Set up SimulationTime
        SimulationTime::Instance()->SetEndTimeAndNumberOfTimeSteps(0.5, 6);

        // Set up mesh
        MutableMesh<2,2> mesh;
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/disk_522_elements");
        mesh.ConstructFromMeshReader(mesh_reader);

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

        // Set up cell population
        MeshBasedCellPopulation<2> cell_population(mesh, cells);

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

        // Create a single PDE and pass to the handler
        // Note SimplePdeForTesting wouldnt work as theres no solution with Neuman conditions.
        // Also note that when using Neuman conditions the only solution that works is u=0
        CellwiseSourcePde<2> pde(cell_population, 0.0);
        ConstBoundaryCondition<2> bc(0.0);
        PdeAndBoundaryConditions<2> pde_and_bc(&pde, &bc, true);
        pde_and_bc.SetDependentVariableName("variable");

        // For coverage, provide an initial guess for the solution
        std::vector<double> data(mesh.GetNumNodes()+1);
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            data[i] = 1.0;
        }

        Vec vector = PetscTools::CreateVec(data);
        pde_and_bc.SetSolution(vector);

        pde_handler.AddPdeAndBc(&pde_and_bc);

        // Open result file ourselves
        OutputFileHandler output_file_handler("TestWritePdeSolution", false);
        pde_handler.mpVizPdeSolutionResultsFile = output_file_handler.OpenOutputFile("results.vizpdesolution");

        // Solve PDE (set sampling timestep multiple to be large doesn't do anything as always output on 1st timestep)
        pde_handler.SolvePdeAndWriteResultsToFile(10);

        // Close result file ourselves
        pde_handler.mpVizPdeSolutionResultsFile->close();

        // Check the correct solution was obtained
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            // Test that PDE solver is working correctly
            TS_ASSERT_DELTA(cell_iter->GetCellData()->GetItem("variable"), 0.0, 0.02);
        }
    }

    void TestSolvePdeAndWriteResultsToFileAndGetPDESolutionAtPointWithoutCoarsePdeMeshDirichlet() throw(Exception)
    {
        EXIT_IF_PARALLEL;

        // Set up SimulationTime
        SimulationTime::Instance()->SetEndTimeAndNumberOfTimeSteps(0.5, 6);

        // Set up mesh
        MutableMesh<2,2> mesh;
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/disk_522_elements");
        mesh.ConstructFromMeshReader(mesh_reader);

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

        // Set up cell population
        MeshBasedCellPopulation<2> cell_population(mesh, cells);

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

        // Create a single PDE and pass to the handler
        SimplePdeForTesting pde;
        ConstBoundaryCondition<2> bc(1.0);
        PdeAndBoundaryConditions<2> pde_and_bc(&pde, &bc, false);
        pde_and_bc.SetDependentVariableName("variable");

        // For coverage, provide an initial guess for the solution
        std::vector<double> data(mesh.GetNumNodes());
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            data[i] = 1.0;
        }

        Vec vector = PetscTools::CreateVec(data);
        pde_and_bc.SetSolution(vector);

        pde_handler.AddPdeAndBc(&pde_and_bc);

        // Open result file ourselves
        OutputFileHandler output_file_handler("TestWritePdeSolution", false);
        pde_handler.mpVizPdeSolutionResultsFile = output_file_handler.OpenOutputFile("results.vizpdesolution");

        // Solve PDE (set sampling timestep multiple to be large doesn't do anything as always output on 1st timestep)
        pde_handler.SolvePdeAndWriteResultsToFile(10);

        // Close result file ourselves
        pde_handler.mpVizPdeSolutionResultsFile->close();

        // Test that this is correct by comparing with an existing results file
        std::string results_dir = output_file_handler.GetOutputDirectoryFullPath();

        NumericFileComparison comparison(results_dir + "results.vizpdesolution", "cell_based/test/data/TestCellBasedPdeHandler/results.vizpdesolution");
        TS_ASSERT(comparison.CompareFiles());

        // Check the correct solution was obtained
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            double radius = norm_2(cell_population.GetLocationOfCellCentre(*cell_iter));
            double analytic_solution = 1.0 - 0.25*(1 - pow(radius,2.0));

            // Test that PDE solver is working correctly
            TS_ASSERT_DELTA(cell_iter->GetCellData()->GetItem("variable"), analytic_solution, 0.02);
        }

        // Now check the GetPdeSolutionAtPoint method

        // First loop over nodes and check it works
        // Check the correct solution was obtained
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            double cell_data_solution(cell_iter->GetCellData()->GetItem("variable"));
            c_vector<double,2> cell_location = cell_population.GetLocationOfCellCentre(*cell_iter);

            TS_ASSERT_DELTA(pde_handler.GetPdeSolutionAtPoint(cell_location,"variable"), cell_data_solution, 1e-6);
        }

        // Now choose some other points

        // Centre
        c_vector<double,2> point;
        point(0) = 0.0;
        point(1) = 0.0;

        TS_ASSERT_DELTA(pde_handler.GetPdeSolutionAtPoint(point,"variable"), 0.75, 0.01);

        // Cover exception
        TS_ASSERT_THROWS_CONTAINS(pde_handler.GetPdeSolutionAtPoint(point, "not_a_var"),
                                  "There is no PDE with that variable.");

        // Random point
        point(0) = 0.5;
        point(1) = 0.5;

//.........这里部分代码省略.........

    void TestRadialSloughingCellKillerMethods() throw(Exception)
    {
        // Create mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_128_elements");
        MutableMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);
        mesh.Translate(-0.5,-0.5);

        // Get centre of mesh (we know it's at the origin, really)
        c_vector<double,2> centre(2);
        centre[0] = 0.0;
        centre[1] = 0.0;
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            centre += mesh.GetNode(i)->rGetLocation();
        }
        centre = centre/mesh.GetNumNodes();

        // Choose radius of cell killer
        double radius = 0.4;

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

        // Create cell population
        MeshBasedCellPopulation<2> cell_population(mesh, cells);

        // Create cell killer and kill cells
        RadialSloughingCellKiller radial_cell_killer(&cell_population, centre, radius);
        radial_cell_killer.CheckAndLabelCellsForApoptosisOrDeath();

        // Check that cells were labelled for death correctly
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            double r = norm_2(cell_population.GetLocationOfCellCentre(*cell_iter) - centre);

            if (r > radius)
            {
                TS_ASSERT_EQUALS(cell_iter->IsDead(), true);
            }
            else
            {
                TS_ASSERT_EQUALS(cell_iter->IsDead(), false);
            }
        }

        // Now get rid of dead cells
        cell_population.RemoveDeadCells();

        // Check that we are correctly left with cells inside the circle of death
        for (AbstractCellPopulation<2>::Iterator cell_iter = cell_population.Begin();
             cell_iter != cell_population.End();
             ++cell_iter)
        {
            double r = norm_2(cell_population.GetLocationOfCellCentre(*cell_iter) - centre);
            TS_ASSERT_LESS_THAN_EQUALS(r, radius);
        }
    }

    void TestHeuristicCleanSimpleMesh() throw(Exception)
    {
        {
            TrianglesMeshReader<1,3> mesh_reader("lung/test/airway_generation/data/test_major_airways_mesh");
            MutableMesh<1,3> mesh;
            mesh.ConstructFromMeshReader(mesh_reader);

            //Assign valid radii
            for(unsigned node_index = 0; node_index < mesh.GetNumNodes(); ++node_index)
            {
                 mesh.GetNode(node_index)->rGetNodeAttributes()[0] = 1.0;
            }

            TS_ASSERT_EQUALS(mesh.GetNumNodes(), 6u);
            TS_ASSERT_EQUALS(mesh.GetNumElements(), 5u);

            MajorAirwaysCentreLinesCleaner cleaner(mesh, 0u);
            cleaner.CleanTerminalsHueristic();

            TS_ASSERT_EQUALS(mesh.GetNumNodes(), 6u);
            TS_ASSERT_EQUALS(mesh.GetNumElements(), 5u);

            //Trips an added nodes exception. Not sure if this is needed!
            //NodeMap node_map(mesh.GetNumAllNodes());
            //mesh.ReIndex(node_map);

            TS_ASSERT_DELTA(mesh.GetNode(2)->rGetLocation()[0], 1.6, 1e-6);
            TS_ASSERT_DELTA(mesh.GetNode(3)->rGetLocation()[0], -1.6, 1e-6);
            TS_ASSERT_DELTA(mesh.GetNode(4)->rGetLocation()[1], 1.6, 1e-6);
            TS_ASSERT_DELTA(mesh.GetNode(5)->rGetLocation()[1], -1.6, 1e-6);


            //Uncomment to visualise
//            VtkMeshWriter<1,3> mesh_writer("TestMajorAirwaysCentreLinesCleaner", "heuristic_shorten", false);
//            mesh_writer.WriteFilesUsingMesh(mesh);
        }

        {
            TrianglesMeshReader<1,3> mesh_reader("lung/test/airway_generation/data/test_major_airways_mesh_short_terminals");
            MutableMesh<1,3> mesh;
            mesh.ConstructFromMeshReader(mesh_reader);

            //Assign valid radii
            for(unsigned node_index = 0; node_index < mesh.GetNumNodes(); ++node_index)
            {
                 mesh.GetNode(node_index)->rGetNodeAttributes()[0] = 1.0;
            }

            TS_ASSERT_EQUALS(mesh.GetNumNodes(), 6u);
            TS_ASSERT_EQUALS(mesh.GetNumElements(), 5u);

            MajorAirwaysCentreLinesCleaner cleaner(mesh, 0u);
            cleaner.CleanTerminalsHueristic();

            TS_ASSERT_EQUALS(mesh.GetNumNodes(), 6u);
            TS_ASSERT_EQUALS(mesh.GetNumElements(), 5u);

            //Trips an added nodes exception. Not sure if this is needed!
            //NodeMap node_map(mesh.GetNumAllNodes());
            //mesh.ReIndex(node_map);

            TS_ASSERT_DELTA(mesh.GetNode(2)->rGetLocation()[0], 1.4, 1e-6);
            TS_ASSERT_DELTA(mesh.GetNode(3)->rGetLocation()[0], -1.4, 1e-6);
            TS_ASSERT_DELTA(mesh.GetNode(4)->rGetLocation()[1], 1.4, 1e-6);
            TS_ASSERT_DELTA(mesh.GetNode(5)->rGetLocation()[1], -1.4, 1e-6);

            //Uncomment to visualise
//            VtkMeshWriter<1,3> mesh_writer("TestMajorAirwaysCentreLinesCleaner", "heuristic_lengthen", false);
//            mesh_writer.WriteFilesUsingMesh(mesh);
        }
    }

    void TestCryptProjectionForceWithArchiving() throw (Exception)
    {
        EXIT_IF_PARALLEL;    // Cell-based archiving doesn't work in parallel.

        OutputFileHandler handler("archive", false);    // don't erase contents of folder
        std::string archive_filename = handler.GetOutputDirectoryFullPath() + "crypt_projection_spring_system.arch";

        {
            TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_2_elements");

            MutableMesh<2,2> mesh;
            mesh.ConstructFromMeshReader(mesh_reader);

            SimulationTime::Instance()->SetEndTimeAndNumberOfTimeSteps(1.0,1);

            std::vector<CellPtr> cells;
            boost::shared_ptr<AbstractCellMutationState> p_state(new WildTypeCellMutationState);
            boost::shared_ptr<AbstractCellProperty> p_stem_type(new StemCellProliferativeType);

            for (unsigned i=0; i<mesh.GetNumNodes(); i++)
            {
                FixedDurationGenerationBasedCellCycleModel* p_model = new FixedDurationGenerationBasedCellCycleModel();
                CellPtr p_cell(new Cell(p_state, p_model));
                p_cell->SetCellProliferativeType(p_stem_type);
                p_cell->SetBirthTime(-50.0);
                cells.push_back(p_cell);
            }

            MeshBasedCellPopulation<2> crypt(mesh, cells);
            WntConcentration<2>::Instance()->SetCryptProjectionParameterA(1.0);
            WntConcentration<2>::Instance()->SetCryptProjectionParameterB(2.0);

            // Create force object
            CryptProjectionForce crypt_projection_force;

            TS_ASSERT_DELTA(crypt_projection_force.GetWntChemotaxisStrength(), 100.0, 1e-6);
            crypt_projection_force.SetWntChemotaxisStrength(15.0);

            std::ofstream ofs(archive_filename.c_str());
            boost::archive::text_oarchive output_arch(ofs);

            // Serialize via pointer
            CryptProjectionForce* const p_crypt_projection_force = &crypt_projection_force;

            p_crypt_projection_force->SetCutOffLength(1.1);

            output_arch << p_crypt_projection_force;
            WntConcentration<2>::Destroy();
        }

        {
            ArchiveLocationInfo::SetMeshPathname("mesh/test/data/", "square_2_elements");

            // Create an input archive
            std::ifstream ifs(archive_filename.c_str(), std::ios::binary);
            boost::archive::text_iarchive input_arch(ifs);

            CryptProjectionForce* p_crypt_projection_force;

            // Restore from the archive
            input_arch >> p_crypt_projection_force;

            // Test the member data
            TS_ASSERT_EQUALS(p_crypt_projection_force->mUseCutOffLength, true);
            TS_ASSERT_DELTA(p_crypt_projection_force->GetA(), 1.0, 1e-12);
            TS_ASSERT_DELTA(p_crypt_projection_force->GetB(), 2.0, 1e-12);
            TS_ASSERT_DELTA(p_crypt_projection_force->GetWntChemotaxisStrength(), 15.0, 1e-6);

            delete p_crypt_projection_force;
        }
    }

    void TestCellwiseDataGradientFineMesh() throw(Exception)
    {
        // Create a mesh: [0,2]x[0,2]
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_4096_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);

        //////////////////////////////////
        // C(x,y) = const
        //////////////////////////////////
        cell_population.SetDataOnAllCells("const", 1.0);

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

        // Check gradient
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            TS_ASSERT_DELTA(gradient.rGetGradient(i)(0), 0.0, 1e-9);
            TS_ASSERT_DELTA(gradient.rGetGradient(i)(1), 0.0, 1e-9);
        }

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

        // Check gradient
        gradient.SetupGradients(cell_population, "x-y");
        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), -1.0, 1e-9);
        }

        // Check gradient - here there is some numerical error
        gradient.SetupGradients(cell_population, "x^2 - y^2");
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            double x = mesh.GetNode(i)->rGetLocation()[0];
            double y = mesh.GetNode(i)->rGetLocation()[1];

            double tol = 0.3;
            if (x==0 || x==2 || y==0 || y==2) //ie on boundary
            {
                tol = 0.6;
            }

            TS_ASSERT_DELTA(gradient.rGetGradient(i)(0),  2*x, tol);
            TS_ASSERT_DELTA(gradient.rGetGradient(i)(1), -2*y, tol);
        }
    }