C++ TS_ASSERT_DELTA函数代码示例

    void TestCellwiseSourcePdeArchiving() throw(Exception)
    {
        EXIT_IF_PARALLEL;

        // Set up simulation time
        SimulationTime::Instance()->SetEndTimeAndNumberOfTimeSteps(1.0, 1);

        // Set up cell population
        HoneycombMeshGenerator generator(5, 5, 0);
        MutableMesh<2,2>* p_mesh = generator.GetMesh();
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumNodes());
        MeshBasedCellPopulation<2> cell_population(*p_mesh, cells);

        FileFinder archive_dir("archive", RelativeTo::ChasteTestOutput);
        std::string archive_file = "CellwiseSourcePde.arch";
        ArchiveLocationInfo::SetMeshFilename("CellwiseSourcePde");

        {
            // Create a PDE object
            AbstractLinearEllipticPde<2,2>* const p_pde = new CellwiseSourcePde<2>(cell_population, 0.05);

            // Create output archive and archive PDE object
            ArchiveOpener<boost::archive::text_oarchive, std::ofstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_oarchive* p_arch = arch_opener.GetCommonArchive();
            (*p_arch) << p_pde;

            delete p_pde;
        }

        {
            AbstractLinearEllipticPde<2,2>* p_pde;

            // Create an input archive and restore PDE object from archive
            ArchiveOpener<boost::archive::text_iarchive, std::ifstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_iarchive* p_arch = arch_opener.GetCommonArchive();
            (*p_arch) >> p_pde;

            // Test that the PDE and its member variables were archived correctly
            TS_ASSERT(dynamic_cast<CellwiseSourcePde<2>*>(p_pde) != NULL);

            CellwiseSourcePde<2>* p_static_cast_pde = static_cast<CellwiseSourcePde<2>*>(p_pde);
            TS_ASSERT_DELTA(p_static_cast_pde->GetCoefficient(), 0.05, 1e-6);
            TS_ASSERT_EQUALS(p_static_cast_pde->mrCellPopulation.GetNumRealCells(), 25u);

            // Avoid memory leaks
            delete &(p_static_cast_pde->mrCellPopulation);
            delete p_pde;
        }
    }

    void TestArchiveAdhesionPottsUpdateRule() throw(Exception)
    {
        OutputFileHandler handler("archive", false);
        std::string archive_filename = handler.GetOutputDirectoryFullPath() + "AdhesionPottsUpdateRule.arch";

        {
            AdhesionPottsUpdateRule<2> update_rule;

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

            // Set member variables
            update_rule.SetCellCellAdhesionEnergyParameter(0.5);
            update_rule.SetCellBoundaryAdhesionEnergyParameter(0.6);

            // Serialize via pointer to most abstract class possible
            AbstractPottsUpdateRule<2>* const p_update_rule = &update_rule;
            output_arch << p_update_rule;
        }

        {
            AbstractPottsUpdateRule<2>* p_update_rule;

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

            // Restore from the archive
            input_arch >> p_update_rule;

            // Test the member data
            TS_ASSERT_DELTA((static_cast<AdhesionPottsUpdateRule<2>*>(p_update_rule))->GetCellCellAdhesionEnergyParameter(), 0.5, 1e-6);
            TS_ASSERT_DELTA((static_cast<AdhesionPottsUpdateRule<2>*>(p_update_rule))->GetCellBoundaryAdhesionEnergyParameter(), 0.6, 1e-6);

            // Tidy up
            delete p_update_rule;
        }
    }

    void TestDerivedQuantities() throw (Exception)
    {
        ParameterisedOde ode;
        boost::shared_ptr<const AbstractOdeSystemInformation> p_info = ode.GetSystemInformation();

        TS_ASSERT_EQUALS(ode.GetNumberOfDerivedQuantities(), 1u);
        TS_ASSERT_EQUALS(ode.HasDerivedQuantity("2a_plus_y"), true);
        TS_ASSERT_EQUALS(ode.HasDerivedQuantity("Not_there"), false);
        TS_ASSERT_EQUALS(ode.GetDerivedQuantityIndex("2a_plus_y"), 0u);
        TS_ASSERT_EQUALS(ode.GetDerivedQuantityUnits(0u), "dimensionless");
        TS_ASSERT_EQUALS(ode.rGetDerivedQuantityNames().size(), 1u);
        TS_ASSERT_EQUALS(ode.rGetDerivedQuantityUnits().size(), 1u);
        TS_ASSERT_EQUALS(ode.rGetDerivedQuantityNames()[0], "2a_plus_y");
        TS_ASSERT_EQUALS(ode.rGetDerivedQuantityUnits()[0], "dimensionless");

        TS_ASSERT_EQUALS(p_info->HasDerivedQuantity("2a_plus_y"), true);
        TS_ASSERT_EQUALS(p_info->HasDerivedQuantity("Not_there"), false);
        TS_ASSERT_EQUALS(p_info->GetDerivedQuantityIndex("2a_plus_y"), 0u);
        TS_ASSERT_EQUALS(p_info->GetDerivedQuantityUnits(0u), "dimensionless");
        TS_ASSERT_EQUALS(p_info->rGetDerivedQuantityNames().size(), 1u);
        TS_ASSERT_EQUALS(p_info->rGetDerivedQuantityUnits().size(), 1u);
        TS_ASSERT_EQUALS(p_info->rGetDerivedQuantityNames()[0], "2a_plus_y");
        TS_ASSERT_EQUALS(p_info->rGetDerivedQuantityUnits()[0], "dimensionless");

        TS_ASSERT_EQUALS(ode.HasAnyVariable("2a_plus_y"), true);
        TS_ASSERT_EQUALS(ode.HasAnyVariable("Not_there"), false);
        TS_ASSERT_EQUALS(ode.GetAnyVariableIndex("2a_plus_y"), 2u);
        TS_ASSERT_EQUALS(ode.GetAnyVariableUnits(2u), "dimensionless");
        TS_ASSERT_EQUALS(p_info->GetAnyVariableIndex("2a_plus_y"), 2u);
        TS_ASSERT_EQUALS(p_info->GetAnyVariableUnits(2u), "dimensionless");

        std::vector<double> derived = ode.ComputeDerivedQuantitiesFromCurrentState(0.0);
        double a = ode.GetParameter(0);
        TS_ASSERT_EQUALS(a, 0.0);
        TS_ASSERT_DELTA(derived[0], 2*a, 1e-4);
        TS_ASSERT_DELTA(ode.GetAnyVariable(2u, 0.0), 2*a, 1e-4);
        TS_ASSERT_DELTA(ode.GetAnyVariable(2u, 0.0, &derived), 2*a, 1e-4);
        a = 1.0;
        ode.SetParameter(0, a);
        derived = ode.ComputeDerivedQuantities(0.0, ode.GetInitialConditions());
        TS_ASSERT_DELTA(derived[0], 2*a, 1e-4);
        double y = 10.0;
        ode.SetStateVariable(0, y);
        derived = ode.ComputeDerivedQuantitiesFromCurrentState(0.0);
        TS_ASSERT_DELTA(derived[0], 2*a+y, 1e-4);
        TS_ASSERT_DELTA(ode.GetAnyVariable(2u, 1.0/* ignored for this ODE */), 2*a+y, 1e-4);

        // Exceptions
        TS_ASSERT_THROWS_THIS(ode.GetDerivedQuantityIndex("Missing"), "No derived quantity named 'Missing'.");
        TS_ASSERT_THROWS_THIS(ode.GetDerivedQuantityUnits(1u), "The index passed in must be less than the number of derived quantities.");

        TwoDimOdeSystem ode2;
        TS_ASSERT_THROWS_THIS(ode2.ComputeDerivedQuantitiesFromCurrentState(0.0),
                              "This ODE system does not define derived quantities.");
        std::vector<double> doesnt_matter;
        TS_ASSERT_THROWS_THIS(ode2.ComputeDerivedQuantities(0.0, doesnt_matter),
                              "This ODE system does not define derived quantities.");
    }

    void TestMirams2010WntOdeSystemSetup() throw(Exception)
    {
#ifdef CHASTE_CVODE
        double wnt_level = 0.5;
        boost::shared_ptr<AbstractCellMutationState> p_state(new WildTypeCellMutationState);
        Mirams2010WntOdeSystem wnt_system(wnt_level, p_state);
        // Solve system using CVODE solver
        // Matlab's strictest bit uses 0.01 below and relaxes it on flatter bits.
        double h_value = 0.1;

        CvodeAdaptor cvode_solver;

        OdeSolution solutions;
        //OdeSolution solutions2;

        std::vector<double> initial_conditions = wnt_system.GetInitialConditions();

        double start_time, end_time, elapsed_time = 0.0;
        start_time = (double) std::clock();
        solutions = cvode_solver.Solve(&wnt_system, initial_conditions, 0.0, 100.0, h_value, h_value);
        end_time = (double) std::clock();
        elapsed_time = (end_time - start_time)/(CLOCKS_PER_SEC);
        std::cout << "1. Cvode Elapsed time = " << elapsed_time << " secs for 100 hours\n";

        // Test solutions are OK for a small time increase...
        int end = solutions.rGetSolutions().size() - 1;
        // Tests the simulation is ending at the right time...(going into S phase at 7.8 hours)
        TS_ASSERT_DELTA(solutions.rGetTimes()[end], 100, 1e-2);

        // Decent results
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][0], 67.5011, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][1], 67.5011, 1e-4);
        TS_ASSERT_DELTA(solutions.rGetSolutions()[end][2], wnt_level, 1e-4);
#else
        std::cout << "CVODE is not enabled. " << std::endl;
        std::cout << "If required please install and alter your hostconfig settings to switch on chaste support." << std::endl;
#endif //CHASTE_CVODE
    }

    void TestCellwiseSourcePdeMethods() throw(Exception)
    {
        EXIT_IF_PARALLEL;

        // Set up cell population
        HoneycombMeshGenerator generator(5, 5, 0);
        MutableMesh<2,2>* p_mesh = generator.GetMesh();
        std::vector<CellPtr> cells;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasic(cells, p_mesh->GetNumNodes());
        MeshBasedCellPopulation<2> cell_population(*p_mesh, cells);

        // Create a PDE object
        CellwiseSourcePde<2> pde(cell_population, 0.05);

        // Test that the member variables have been initialised correctly
        TS_ASSERT_EQUALS(&(pde.rGetCellPopulation()), &cell_population);
        TS_ASSERT_DELTA(pde.GetCoefficient(), 0.05, 1e-6);

        // Test methods
        Node<2>* p_node = cell_population.GetNodeCorrespondingToCell(*(cell_population.Begin()));
        TS_ASSERT_DELTA(pde.ComputeLinearInUCoeffInSourceTermAtNode(*p_node), 0.05, 1e-6);

        ChastePoint<2> point;
        c_matrix<double,2,2> diffusion_matrix = pde.ComputeDiffusionTerm(point);
        for (unsigned i=0; i<2; i++)
        {
            for (unsigned j=0; j<2; j++)
            {
                double value = 0.0;
                if (i == j)
                {
                    value = 1.0;
                }
                TS_ASSERT_DELTA(diffusion_matrix(i,j), value, 1e-6);
            }
        }
    }

    void TestMonodomainTissueGetCardiacCell() throw(Exception)
    {
        HeartConfig::Instance()->Reset();
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(1.0, 1.0); // [0,1] with h=1.0, ie 2 node mesh

        MyCardiacCellFactory cell_factory;
        cell_factory.SetMesh(&mesh);

        MonodomainTissue<1> monodomain_tissue( &cell_factory );

        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(0))
        {
            AbstractCardiacCellInterface* cell = monodomain_tissue.GetCardiacCell(0);
            TS_ASSERT_DELTA(cell->GetStimulus(0.001),-80,1e-10);
        }

        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(1))
        {
            AbstractCardiacCellInterface* cell = monodomain_tissue.GetCardiacCell(1);
            TS_ASSERT_DELTA(cell->GetStimulus(0.001),0,1e-10);
        }
    }

        void test_initial_settings()
        {
            TS_ASSERT_EQUALS(ch_->num_links(), NUM_MODULES);
            for (int i=0; i < ch_->num_links(); i++)
            {
                TS_ASSERT_DELTA(ch_->get_link(i).get_q(), INITIAL_Q, EPS);
            }

            for (int i=0; i < ch_->num_links(); i++)
            {
                ch_->get_link(i).set_q(0.0);
                TS_ASSERT_DELTA(ch_->get_link(i).get_q(), 0.0, EPS);
                TS_ASSERT((ch_->get_link(i).get_R_jts()).isApprox(Eigen::Matrix3d::Identity()));
            }
            /* Angles are now 0, so all of the modules should be
             * at the "initial_rotaiton" rotation (as long as
             * R_jts is identity for all of the modules)
             */
            for (int i=0; i < ch_->num_links(); i++)
            {
                TS_ASSERT(ch_->get_current_R(i).isApprox(ch_->get_link(0).get_init_rotation()));
            }
        } 

    void TestCalculateLobeVolume()
    {
    #if defined(CHASTE_VTK) && ( (VTK_MAJOR_VERSION >= 5 && VTK_MINOR_VERSION >= 6) || VTK_MAJOR_VERSION >= 6)
       EXIT_IF_PARALLEL;

       vtkSmartPointer<vtkPolyData> sphere = CreateSphere(100);

       AirwayGenerator generator(sphere);

       //The sphere is coarsely meshed, hence relatively large tolerance
       TS_ASSERT_DELTA(generator.CalculateLobeVolume(), 4.0/3.0*M_PI, 1e-2);

    #endif
    }

    void TestXaxisRotation3DWithMethod()
    {
        TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_136_elements");
        TetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        ChastePoint<3> corner_before = mesh.GetNode(6)->GetPoint();
        TS_ASSERT_EQUALS(corner_before[0], 1.0);
        TS_ASSERT_EQUALS(corner_before[1], 1.0);
        TS_ASSERT_EQUALS(corner_before[2], 1.0);

        double mesh_volume = mesh.GetVolume();

        mesh.RotateX(M_PI/2.0);

        double new_mesh_volume = mesh.GetVolume();
        TS_ASSERT_DELTA(mesh_volume,new_mesh_volume,1e-6);

        ChastePoint<3> corner_after = mesh.GetNode(6)->GetPoint();
        TS_ASSERT_DELTA(corner_after[0],  1.0, 1e-7);
        TS_ASSERT_DELTA(corner_after[1],  1.0, 1e-7);
        TS_ASSERT_DELTA(corner_after[2], -1.0, 1e-7);
    }

    void TestHeatEquationForCoupledOdeSystem()
    {
        // Create PDE system object
        HeatEquationForCoupledOdeSystem<1> pde;

        ChastePoint<1> x(1.0);

        TS_ASSERT_DELTA(pde.ComputeDuDtCoefficientFunction(x,0), 1.0, 1e-6);

        c_vector<double,1> pde_solution;
        pde_solution(0) = 4.0;

        std::vector<double> ode_solution(1);
        ode_solution[0] = 5.0;

        TS_ASSERT_DELTA(pde.ComputeSourceTerm(x, pde_solution, ode_solution, 0), 0.0, 1e-6);

        Node<1> node(0);
        TS_ASSERT_DELTA(pde.ComputeSourceTermAtNode(node, pde_solution, ode_solution, 0), 0.0, 1e-6);

        c_matrix<double,1,1> diffusion_term = pde.ComputeDiffusionTerm(x, 0);
        TS_ASSERT_DELTA(diffusion_term(0,0), 1.0, 1e-6);
    }

    void TestLoadAbstractOdeSystem()
    {
        TwoDimOdeSystem ode;

        TS_ASSERT_EQUALS( ode.GetNumberOfStateVariables(), 2U );

        // Read archive from previous test
        OutputFileHandler handler("archive", false);
        std::string archive_filename;
        archive_filename = handler.GetOutputDirectoryFullPath() + "ode.arch";

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

        input_arch >> ode;

        TS_ASSERT_EQUALS( ode.GetNumberOfStateVariables(), 2U );

        std::vector<double>& r_state_variables = ode.rGetStateVariables();

        TS_ASSERT_DELTA(r_state_variables[0], 7.0, 1e-12);
        TS_ASSERT_DELTA(r_state_variables[1], 8.0, 1e-12);
    }

    /*
     * To visualize the results, open a new terminal, {{{cd}}} to the Chaste directory,
     * then {{{cd}}} to {{{anim}}}. Then do: {{{java Visualize2dVertexCells /tmp/$USER/testoutput/CellBasedDemo1/results_from_time_0}}}.
     * We may have to do: {{{javac Visualize2dVertexCells.java}}} beforehand to create the
     * java executable.
     *
     * EMPTYLINE
     *
     * The {{{make_a_movie}}} script can be used to generate a video based on the results of your simulation.
     * To do this, first visualize the results using {{{Visualize2dVertexCells}}} as described above. Click
     * on the box marked "Output" and play through the whole simulation to generate a sequence of {{{.png}}}
     * images, one for each time step. Next, still in the {{{anim}}} folder, do: {{{./make_a_movie}}}.
     * This reads in the {{{.png}}} files and creates a video file called {{{simulation.mpeg}}}.
     *
     * EMPTYLINE
     *
     * Results can also be visualized using Paraview. See the UserTutorials/VisualizingWithParaview tutorial for more information.
     *
     * EMPTYLINE
     *
     * == Test 2 - basic node-based simulation ==
     *
     * We next show how to modify the previous test to implement a 'node-based' simulation,
     * in which cells are represented by overlapping spheres (actually circles, since we're
     * in 2D).
     */
    void TestNodeBasedMonolayer() throw (Exception)
    {
        /* We now need to create a {{{NodesOnlyMesh}}} we do this by first creating a {{{MutableMesh}}}
         * and passing this to a helper method {{{ConstructNodesWithoutMesh}}} along with a interaction cut off length
         * that defines the connectivity in the mesh.
         */
        HoneycombMeshGenerator generator(2, 2); //**Changed**//
        MutableMesh<2,2>* p_generating_mesh = generator.GetMesh(); //**Changed**//
        NodesOnlyMesh<2> mesh; //**Changed**//
        mesh.ConstructNodesWithoutMesh(*p_generating_mesh, 1.5); //**Changed**//

        /* We create the cells as before, only this time we need one cell per node.*/
        std::vector<CellPtr> cells;
        MAKE_PTR(TransitCellProliferativeType, p_transit_type);
        CellsGenerator<StochasticDurationCellCycleModel, 2> cells_generator;
        cells_generator.GenerateBasicRandom(cells, mesh.GetNumNodes(), p_transit_type); //**Changed**//

        /* This time we create a {{{NodeBasedCellPopulation}}} as we are using a {{{NodesOnlyMesh}}}.*/
        NodeBasedCellPopulation<2> cell_population(mesh, cells);//**Changed**//

        /* We create an {{{OffLatticeSimulation}}} object as before, all we change is the output directory
         * and output results more often as a larger default timestep is used for these simulations. */
        OffLatticeSimulation<2> simulator(cell_population);
        simulator.SetOutputDirectory("CellBasedDemo2"); //**Changed**//
        simulator.SetSamplingTimestepMultiple(12); //**Changed**//
        simulator.SetEndTime(20.0);

        /* We use a different {{{Force}}} which is suitable for node based simulations.
         */
        MAKE_PTR(RepulsionForce<2>, p_force); //**Changed**//
        simulator.AddForce(p_force);

        /* In all types of simulation you may specify how cells are removed from the simulation by specifying
         * a {{{CellKiller}}}. You create these in the same was as the {{{Force}}} and pass them to the {{{CellBasedSimulation}}}.
         * Note that here the constructor for {{{RandomCellKiller}}} requires some arguments to be passed to it, therefore we use the
         * {{{MAKE_PTR_ARGS}}} macro.
         */
        MAKE_PTR_ARGS(RandomCellKiller<2>, p_cell_killer, (&cell_population, 0.01)); //**Changed**//
        simulator.AddCellKiller(p_cell_killer);

        /* Again we call the {{{Solve}}} method on the simulation to run the simulation.*/
        simulator.Solve();

        /* The next two lines are for test purposes only and are not part of this tutorial.
         * Again, we are checking that we reached the end time of the simulation
         * with the correct number of cells.
         */
        TS_ASSERT_EQUALS(cell_population.GetNumRealCells(), 7u);
        TS_ASSERT_DELTA(SimulationTime::Instance()->GetTime(), 20.0, 1e-10);
    }

    void TestArchiving() throw (Exception)
    {
        EXIT_IF_PARALLEL;

        OutputFileHandler handler("TestElementAttributes", false);
        std::string archive_filename = handler.GetOutputDirectoryFullPath() + "element_attributes.arch";

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

            ElementAttributes<2,2>* p_element_attributes = new ElementAttributes<2,2>();

            p_element_attributes->AddAttribute(2.34);
            p_element_attributes->AddAttribute(5.67);

            ElementAttributes<2,2>* const p_const_element_attributes = p_element_attributes;

            output_arch << p_const_element_attributes;

            delete p_element_attributes;
        }

        {
            std::ifstream ifs(archive_filename.c_str(), std::ios::binary);
            boost::archive::text_iarchive input_arch(ifs);

            ElementAttributes<2,2>* p_element_attributes;
            input_arch >> p_element_attributes;

            TS_ASSERT_EQUALS(p_element_attributes->rGetAttributes().size(), 2u);
            TS_ASSERT_DELTA(p_element_attributes->rGetAttributes()[0], 2.34, 1e-10);
            TS_ASSERT_DELTA(p_element_attributes->rGetAttributes()[1], 5.67, 1e-10);

            delete p_element_attributes;
        }
    }

    // Testing Load() (based on previous two tests)
    void TestLoad() throw (Exception)
    {
        EXIT_IF_PARALLEL;    // Cell based archiving doesn't work in parallel.

        // Load the simulation from the TestSave method above and
        // run it from 0.1 to 1.0
        OffLatticeSimulation<2>* p_simulator1;
        p_simulator1 = CellBasedSimulationArchiver<2, OffLatticeSimulation<2> >::Load("TestOffLatticeSimulationWithNodeBasedCellPopulationSaveAndLoad", 0.1);

        // Test that the numerical method was archived correctly
        boost::shared_ptr<AbstractNumericalMethod<2, 2> > p_method = p_simulator1->GetNumericalMethod();
        TS_ASSERT_EQUALS(p_method->HasAdaptiveTimestep(), true);

        p_simulator1->SetEndTime(1.0);
        p_simulator1->Solve();

        // Save, then reload and run from 1.0 to 2.5
        CellBasedSimulationArchiver<2, OffLatticeSimulation<2> >::Save(p_simulator1);
        OffLatticeSimulation<2>* p_simulator2
            = CellBasedSimulationArchiver<2, OffLatticeSimulation<2> >::Load("TestOffLatticeSimulationWithNodeBasedCellPopulationSaveAndLoad", 1.0);

        p_simulator2->SetEndTime(2.5);
        p_simulator2->Solve();

        // These results are from time 2.5 in TestStandardResultForArchivingTestBelow() (above!)
        std::vector<double> node_3_location = p_simulator2->GetNodeLocation(3);
        TS_ASSERT_DELTA(node_3_location[0], mNode3x, 1e-4);
        TS_ASSERT_DELTA(node_3_location[1], mNode3y, 1e-4);

        std::vector<double> node_4_location = p_simulator2->GetNodeLocation(4);
        TS_ASSERT_DELTA(node_4_location[0], mNode4x, 1e-4);
        TS_ASSERT_DELTA(node_4_location[1], mNode4y, 1e-4);

        // Tidy up
        delete p_simulator1;
        delete p_simulator2;
    }

    void TestSimpleSystemWithOrderSwapped()
    {
        // The solution should be the same, but we'll have to construct a new CombinedOdeSystemInformation
        // object, because the order of subsystems has changed.

        SimpleOde1 ode_for_y; // dy/dt = x
        SimpleOde2 ode_for_x; // dx/dt = -y

        std::vector<AbstractOdeSystem*> ode_systems;
        ode_systems.push_back(&ode_for_x);
        ode_systems.push_back(&ode_for_y);

        // Create combined ODE system
        CombinedOdeSystem combined_ode_system(ode_systems);

        // Tell the combined ODE system which state variables in the first ODE system
        // correspond to which parameters in the second ODE system...
        std::map<unsigned, unsigned> variable_parameter_map;
        variable_parameter_map[0] = 0;

        combined_ode_system.Configure(variable_parameter_map, &ode_for_y, &ode_for_x);

        // ...and vice versa (we can re-use the map in this case)
        combined_ode_system.Configure(variable_parameter_map, &ode_for_x, &ode_for_y);

        // Test solving the combined system.
        // This is dy/dt = x, dx/dt = -y, y(0) = 0, x(0) = 1.
        // The analytic solution is y = sin(t), x = cos(t).
        EulerIvpOdeSolver solver;
        OdeSolution solutions;
        double h = 0.01;
        std::vector<double> inits = combined_ode_system.GetInitialConditions();
        solutions = solver.Solve(&combined_ode_system, inits, 0.0, 2.0, h, h);
        double global_error = 0.5*(exp(2.0)-1)*h;
        TS_ASSERT_DELTA(solutions.rGetSolutions().back()[1], sin(2.0), global_error);
        TS_ASSERT_DELTA(solutions.rGetSolutions().back()[0], cos(2.0), global_error);
    }