C++ MutableMesh类代码示例

    void TestConstructorWithCellPopulation1d() throw (Exception)
    {
        // Create 1D cell population
        unsigned num_cells = 5;
        MutableMesh<1,1> mesh;
        mesh.ConstructLinearMesh(num_cells-1);

        std::vector<CellPtr> cells1d;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator1d;
        cells_generator1d.GenerateBasic(cells1d, mesh.GetNumNodes());
        MeshBasedCellPopulation<1> crypt1d(mesh, cells1d);

        // Pass this cell population to a boundary condition object
        CryptSimulationBoundaryCondition<1> boundary_condition1d(&crypt1d);

        // Test access to the cell population
        const AbstractCellPopulation<1>* p_population1d = boundary_condition1d.GetCellPopulation();
        TS_ASSERT(p_population1d != NULL);
    }

    void TestSloughingCellKillerIn3d() throw(Exception)
    {
        // Create 3D mesh
        MutableMesh<3,3> mesh;
        mesh.ConstructCuboid(4, 5, 6);

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

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

        // Create cell killer
        SloughingCellKiller<3> sloughing_cell_killer(&cell_population, 1.0); // number is irrelevent as long as its positive.

        // Check that an exception is thrown, as this method is not yet implemented in 3D
        TS_ASSERT_THROWS_THIS(sloughing_cell_killer.CheckAndLabelCellsForApoptosisOrDeath(), "SloughingCellKiller is not yet implemented in 3D");
    }

    void TestVerifyBoundaryCondition1d() throw(Exception)
    {
        // Create 1D cell population
        unsigned num_cells = 5;
        MutableMesh<1,1> mesh;
        mesh.ConstructLinearMesh(num_cells-1);

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

        // Set up a WntConcentration singleton object
        WntConcentration<1>* p_wnt = WntConcentration<1>::Instance();
        WntConcentration<1>::Instance()->SetType(LINEAR);
        WntConcentration<1>::Instance()->SetCellPopulation(crypt);
        WntConcentration<1>::Instance()->SetCryptLength(crypt.GetWidth(0));
        TS_ASSERT_EQUALS(p_wnt->IsWntSetUp(), true);

        // Store the node locations
        std::map<Node<1>*, c_vector<double, 1> > node_locations_before;
        for (unsigned node_index=0; node_index<mesh.GetNumNodes(); node_index++)
        {
            node_locations_before[crypt.GetNode(node_index)] = crypt.GetNode(node_index)->rGetLocation();
        }

        crypt.GetNode(0)->rGetModifiableLocation()[0] = -0.1;

        // Create a boundary condition object
        CryptSimulationBoundaryCondition<1> boundary_condition(&crypt);

        // Before imposing the boundary condition, it should not be satisfied
        TS_ASSERT_EQUALS(boundary_condition.VerifyBoundaryCondition(), false);

        boundary_condition.ImposeBoundaryCondition(node_locations_before);

        // After imposing the boundary condition, it should be satisfied
        TS_ASSERT_EQUALS(boundary_condition.VerifyBoundaryCondition(), true);
    }

    void TestImposeBoundaryConditionWithNoWnt1d() throw(Exception)
    {
        // Create 1D cell population
        unsigned num_cells = 5;
        MutableMesh<1,1> mesh;
        mesh.ConstructLinearMesh(num_cells-1);

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

        // Store the node locations
        std::map<Node<1>*, c_vector<double, 1> > node_locations_before;
        for (unsigned node_index=0; node_index<mesh.GetNumNodes(); node_index++)
        {
            node_locations_before[crypt.GetNode(node_index)] = crypt.GetNode(node_index)->rGetLocation();
        }

        // Now move the first cell (which should be on x=0, and a stem cell) to the left a bit
        AbstractCellPopulation<1>::Iterator cell_iter = crypt.Begin();

        // Check is a stem cell
        TS_ASSERT_EQUALS(cell_iter->GetCellProliferativeType()->IsType<StemCellProliferativeType>(), true);

        // Check initially at x=0
        TS_ASSERT_DELTA(crypt.GetLocationOfCellCentre(*cell_iter)[0], 0.0, 1e-6);

        // Now move to the left a bit
        crypt.GetNode(0)->rGetModifiableLocation()[0] = -0.1;
        TS_ASSERT_LESS_THAN(crypt.GetLocationOfCellCentre(*cell_iter)[0], 0.0);

        // Pass this cell population to a boundary condition object
        CryptSimulationBoundaryCondition<1> boundary_condition(&crypt);

        // Impose the boundary condition without a Wnt stimulus
        boundary_condition.ImposeBoundaryCondition(node_locations_before);

        /*
         * The first cell should have been pulled back to x=0 since it is a stem cell and
         * there is no Wnt stimulus. It should be unaffected by jiggling, which is not imposed
         * in 1D.
         */
        TS_ASSERT_DELTA(0.0, crypt.GetLocationOfCellCentre(*cell_iter)[0], 1e-3);

        // The nodes should all now be at their original locations
        std::map<Node<1>*, c_vector<double, 1> > node_locations_after;
        for (unsigned node_index=0; node_index<mesh.GetNumNodes(); node_index++)
        {
            node_locations_after[crypt.GetNode(node_index)] = crypt.GetNode(node_index)->rGetLocation();
        }

        for (unsigned node_index=0; node_index<mesh.GetNumNodes(); node_index++)
        {
            TS_ASSERT_DELTA(node_locations_before[crypt.GetNode(node_index)](0), node_locations_after[crypt.GetNode(node_index)](0), 1e-3);
        }
    }

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