void TestTranslation2DWithUblas()
    {
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/2D_0_to_1mm_200_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        double volume = mesh.GetVolume();
        double surface_area = mesh.GetSurfaceArea();
        Node<2>* p_node1 = mesh.GetNode(36);
        ChastePoint<2> point1 = p_node1->GetPoint();
        Node<2>* p_node2 = mesh.GetNode(23);
        ChastePoint<2> point2 = p_node2->GetPoint();

        c_vector<double, 2> old_location1 = point1.rGetLocation();
        c_vector<double, 2> old_location2 = point2.rGetLocation();

        // Set translation Vector
        c_vector<double, 2> trans_vec;
        trans_vec(0) = 2.0;
        trans_vec(1) = 2.0;

        // Translate
        mesh.Translate(trans_vec);
        c_vector<double, 2> new_location1 = point1.rGetLocation();
        c_vector<double, 2> new_location2 = point2.rGetLocation();

        // Check Volume and Surface Area are invariant
        TS_ASSERT_DELTA(mesh.GetVolume(), volume, 1e-6);
        TS_ASSERT_DELTA(mesh.GetSurfaceArea(), surface_area, 1e-6);

        // Spot check a couple of nodes
        TS_ASSERT_DELTA(inner_prod(new_location1-old_location1, trans_vec), 0, 1e-6);
        TS_ASSERT_DELTA(inner_prod(new_location2-old_location2, trans_vec), 0, 1e-6);
    }

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

        TS_ASSERT_DELTA(mesh.GetVolume(), 1.0, 1e-6);
        TS_ASSERT_DELTA(mesh.GetSurfaceArea(), 6.0, 1e-6);

        // Change coordinates

        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            Node<3>* p_node = mesh.GetNode(i);
            ChastePoint<3> point = p_node->GetPoint();
            point.SetCoordinate(0, point[0]*2.0);
            point.SetCoordinate(1, point[1]*2.0);
            point.SetCoordinate(2, point[2]*2.0);
            p_node->SetPoint(point);
        }

        mesh.RefreshMesh();

        TS_ASSERT_DELTA(mesh.GetVolume(), 8.0, 1e-6);
        TS_ASSERT_DELTA(mesh.GetSurfaceArea(), 24.0, 1e-6);
    }

    void TestBathIntracellularStimulation() throw (Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(10.0);  //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainBath1d");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain_bath_1d");

        c_vector<double,1> centre;
        centre(0) = 0.5;
        BathCellFactory<1> cell_factory(-1e6, centre); // stimulates x=0.5 node

        BidomainWithBathProblem<1> bidomain_problem( &cell_factory );

        TrianglesMeshReader<1,1> reader("mesh/test/data/1D_0_to_1_100_elements");
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructFromMeshReader(reader);

        // set the x<0.25 and x>0.75 regions as the bath region
        for(unsigned i=0; i<mesh.GetNumElements(); i++)
        {
            double x = mesh.GetElement(i)->CalculateCentroid()[0];
            if( (x<0.25) || (x>0.75) )
            {
                mesh.GetElement(i)->SetAttribute(HeartRegionCode::GetValidBathId());
            }
        }

        bidomain_problem.SetMesh(&mesh);
        bidomain_problem.Initialise();

        bidomain_problem.Solve();

        Vec sol = bidomain_problem.GetSolution();
        ReplicatableVector sol_repl(sol);

        // test V = 0 for all bath nodes
        for(unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            if(HeartRegionCode::IsRegionBath( mesh.GetNode(i)->GetRegion() )) // bath
            {
                TS_ASSERT_DELTA(sol_repl[2*i], 0.0, 1e-12);
            }
        }

        // test symmetry of V and phi_e
        for(unsigned i=0; i<=(mesh.GetNumNodes()-1)/2; i++)
        {
            unsigned opposite = mesh.GetNumNodes()-i-1;
            TS_ASSERT_DELTA(sol_repl[2*i], sol_repl[2*opposite], 2e-3);      // V
            TS_ASSERT_DELTA(sol_repl[2*i+1], sol_repl[2*opposite+1], 2e-3);  // phi_e
        }

        // a couple of hardcoded values
        TS_ASSERT_DELTA(sol_repl[2*50], 3.7684, 1e-3);
        TS_ASSERT_DELTA(sol_repl[2*70], 5.1777, 1e-3);
    }

    void TestCoverage3d()
    {
        HeartConfig::Instance()->SetSimulationDuration(0.1); //ms
        HeartConfig::Instance()->SetUseStateVariableInterpolation(true);
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.005, 0.01, 0.1);

        TetrahedralMesh<3,3> mesh;
        mesh.ConstructRegularSlabMesh(0.02, 0.02, 0.02, 0.02);

        ZeroStimulusCellFactory<CellLuoRudy1991FromCellML,3> cell_factory;
        MonodomainProblem<3> monodomain_problem( &cell_factory );
        monodomain_problem.SetMesh(&mesh);
        monodomain_problem.Initialise();
        monodomain_problem.Solve();
    }

    void TestRemeshSingleBranch() throw(Exception)
    {
        //Load a single branch mesh file
        TrianglesMeshReader<1,3> reader("mesh/test/data/1D_in_3D_0_to_1mm_10_elements");
        TetrahedralMesh<1,3> mesh;
        mesh.ConstructFromMeshReader(reader);

        //We need to add some attributes to the mesh
        for (TetrahedralMesh<1,3>::NodeIterator iter = mesh.GetNodeIteratorBegin();
             iter != mesh.GetNodeIteratorEnd();
             ++iter)
        {
            iter->AddNodeAttribute(0.05);
        }

        //Create remesher object
        AirwayRemesher remesher(mesh, 0u);

        //Check intermediate elements are removed. Poiseuille resistance of the branch is
        // C*0.1/((0.05^4) = C*1.6 * 10^4

        MutableMesh<1,3> output_mesh_one;
        remesher.Remesh(output_mesh_one, 1e5); //With this tolerance all intermediate nodes should be removed.

        TS_ASSERT_EQUALS(output_mesh_one.GetNumNodes(), 2u);
        TS_ASSERT_EQUALS(output_mesh_one.GetNumElements(), 1u);
        TS_ASSERT_DELTA(output_mesh_one.GetElement(0)->GetAttribute(), 0.05, 1e-6);

        MutableMesh<1,3> output_mesh_two;
        remesher.Remesh(output_mesh_two, 0.8e4); //With this tolerance there should be one intermediate node.

        TS_ASSERT_EQUALS(output_mesh_two.GetNumNodes(), 3u);
        TS_ASSERT_EQUALS(output_mesh_two.GetNumElements(), 2u);
        TS_ASSERT_DELTA(output_mesh_two.GetElement(0)->GetAttribute(), 0.05, 1e-6);

        MutableMesh<1,3> output_mesh_three;
        remesher.Remesh(output_mesh_three, 1.6e3); //With this tolerance there should be ten elements.

        TS_ASSERT_EQUALS(output_mesh_three.GetNumNodes(), 11u);
        TS_ASSERT_EQUALS(output_mesh_three.GetNumElements(), 10u);
        TS_ASSERT_DELTA(output_mesh_three.GetElement(0)->GetAttribute(), 0.05, 1e-6);
        TS_ASSERT_DELTA(output_mesh_three.GetElement(5)->GetAttribute(), 0.05, 1e-6);


        //To visualise
        //VtkMeshWriter<1,3> writer("TestAirwayRemesher", "1D_remeshed");
        //writer.WriteFilesUsingMesh(output_mesh_three);
    }

    void TestElementReactanceAndInertance()
    {
        TetrahedralMesh<1,3> mesh;
        TrianglesMeshReader<1,3> mesh_reader("mesh/test/data/y_branch_3d_mesh");
        mesh.ConstructFromMeshReader(mesh_reader);

        SimpleImpedanceProblem problem(mesh, 0u);

        problem.SetRho(M_PI);
        problem.SetMu(M_PI);

        double l = 2.0;
        double r = 2.0;

        TS_ASSERT_DELTA(problem.CalculateElementResistance(r, l), 8*l/(r*r*r*r), 1e-6);
        TS_ASSERT_DELTA(problem.CalculateElementInertance(r, l), l/(r*r), 1e-6);
    }

    void TestMultipleFrequencies() throw(Exception)
    {
        TetrahedralMesh<1,3> mesh;
        //TrianglesMeshReader<1,3> mesh_reader("mesh/test/data/y_branch_3d_mesh");
        TrianglesMeshReader<1,3> mesh_reader("lung/test/data/TestSubject002");
        mesh.ConstructFromMeshReader(mesh_reader);

        //Scale all radii by 0.7 to give an FRC equivalent lung
        for (TetrahedralMesh<1,3>::NodeIterator node_iter = mesh.GetNodeIteratorBegin();
             node_iter != mesh.GetNodeIteratorEnd();
             ++node_iter)
        {
            node_iter->rGetNodeAttributes()[0] *= 0.7;
        }

        std::vector<double> test_frequencies;
        test_frequencies.push_back(1.0);
        test_frequencies.push_back(2.0);
        test_frequencies.push_back(3.0);
        test_frequencies.push_back(5.0);
        test_frequencies.push_back(10.0);
        test_frequencies.push_back(20.0);
        test_frequencies.push_back(30.0);

        SimpleImpedanceProblem problem(mesh, 0u);
        problem.SetMeshInMilliMetres();
        problem.SetFrequencies(test_frequencies);               //Set & get frequencies for coverage
        std::vector<double>& freqs = problem.rGetFrequencies();

        TS_ASSERT_EQUALS(freqs.size(), 7u);

        problem.Solve();

        std::vector<std::complex<double> > impedances = problem.rGetImpedances();

        TS_ASSERT_EQUALS(impedances.size(), 7u);

        //These are hard coded from previous runs, but are as expected for
        //a patient with moderate to severe asthma
        TS_ASSERT_DELTA(real(impedances[0])*1e-3/98, 8.45, 1e-2);
        TS_ASSERT_DELTA(imag(impedances[0])*1e-3/98, -3.65, 1e-2);
        TS_ASSERT_DELTA(real(impedances[6])*1e-3/98, 5.77, 1e-2);
        TS_ASSERT_DELTA(imag(impedances[6])*1e-3/98, 4.12, 1e-2);
    }

    void RunBenchMark(double h, double dt, double endTime, bool useSvi)
    {
        TetrahedralMesh<3,3> mesh;

        mesh.ConstructRegularSlabMesh(h, 2.0, 0.7, 0.3);

        std::stringstream output_dir;
        output_dir << "Benchmark" << "_h" << h << "_dt" << dt;

        HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");
        HeartConfig::Instance()->SetSimulationDuration(endTime); //ms

        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.005, dt, 0.1);
        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1400); // 1400 1/cm
        HeartConfig::Instance()->SetCapacitance(1); // 1uF/cm^2

        HeartConfig::Instance()->SetVisualizeWithMeshalyzer(false);

        // The Chaste results for the benchmark paper use STATE-VARIABLE INTERPOLATION switched on
        // (see comments above)
        HeartConfig::Instance()->SetUseStateVariableInterpolation(useSvi);

        // Regarding the second paper described above, to run the simulations with ICI, comment out the
        // above line. To run the simulation with operator splitting, or with (full) mass-lumping,
        // comment out the above SVI line and uncomment one of the below. (Note: half-lumping is not
        // available).
        //HeartConfig::Instance()->SetUseMassLumping(true); // what is described as full-lumping in this paper
        //HeartConfig::Instance()->SetUseReactionDiffusionOperatorSplitting(true);

        double long_conductance = 0.17 * 0.62/(0.17+0.62) * 10;    // harmonic mean of 0.17, 0.62 S/m converted to mS/cm
        double trans_conductance = 0.019 * 0.24/(0.019+0.24) * 10; // harmonic mean of 0.019,0.24 S/m converted to mS/cm

        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(long_conductance, trans_conductance, trans_conductance));

        BenchmarkCellFactory cell_factory;

        MonodomainProblem<3> problem( &cell_factory );
        problem.SetMesh(&mesh);

        problem.Initialise();
        problem.SetWriteInfo();
        problem.Solve();
    }

    void TestAcinarImpedance() throw(Exception)
    {
        TetrahedralMesh<1,3> mesh;
        TrianglesMeshReader<1,3> mesh_reader("mesh/test/data/y_branch_3d_mesh");
        mesh.ConstructFromMeshReader(mesh_reader);

        TS_ASSERT_THROWS_CONTAINS(SimpleImpedanceProblem(mesh, 1u), "Outlet node is not a boundary node");

        SimpleImpedanceProblem problem(mesh, 0u);
        problem.rGetMesh(); //for coverage

        unsigned node_index = 3; //Arbitrary terminal node

        problem.SetElastance(2*M_PI/2.0);

        TS_ASSERT_DELTA(real(problem.CalculateAcinusImpedance(mesh.GetNode(node_index), 0.0)), 0.0, 1e-6);
        TS_ASSERT_DELTA(real(problem.CalculateAcinusImpedance(mesh.GetNode(node_index), 1.0)), 0.0, 1e-6);
        TS_ASSERT_DELTA(imag(problem.CalculateAcinusImpedance(mesh.GetNode(node_index), 1.0)), -1.0, 1e-6);
    }

    void TestGetSingleRadiusVector(void) throw(Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/simple_cube");
        TetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        TS_ASSERT_EQUALS(mesh.GetNumElements(),12u);

        PapillaryFibreCalculator calculator(mesh);

        // Call GetRadiusVectors on an element
        unsigned element_index = 0;
        c_vector<double, 3> radius_vector = calculator.GetRadiusVectorForOneElement(element_index);

        // Check they are right
        TS_ASSERT_DELTA(radius_vector[0], -0.275, 1e-9);
        TS_ASSERT_DELTA(radius_vector[1], -0.025, 1e-9);
        TS_ASSERT_DELTA(radius_vector[2], -0.275, 1e-9);
    }

    /**
     * Simple Parabolic PDE u' = del squared u
     *
     * With u = 0 on the boundaries of the unit cube. Subject to the initial
     * condition u(0,x,y,z)=sin( PI x)sin( PI y)sin( PI z).
     */
    void TestSimpleLinearParabolicSolver3DZeroDirich()
    {
        // read mesh on [0,1]x[0,1]x[0,1]
        TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_136_elements");
        TetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        // Instantiate PDE object
        HeatEquation<3> pde;

        // Boundary conditions - zero dirichlet everywhere on boundary
        BoundaryConditionsContainer<3,3,1> bcc;
        bcc.DefineZeroDirichletOnMeshBoundary(&mesh);

        // Solver
        SimpleLinearParabolicSolver<3,3> solver(&mesh,&pde,&bcc);

        /*
         * Choose initial condition sin(x*pi)*sin(y*pi)*sin(z*pi) as
         * this is an eigenfunction of the heat equation.
         */
        std::vector<double> init_cond(mesh.GetNumNodes());
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            double x = mesh.GetNode(i)->GetPoint()[0];
            double y = mesh.GetNode(i)->GetPoint()[1];
            double z = mesh.GetNode(i)->GetPoint()[2];
            init_cond[i] = sin(x*M_PI)*sin(y*M_PI)*sin(z*M_PI);
        }
        Vec initial_condition = PetscTools::CreateVec(init_cond);

        double t_end = 0.1;
        solver.SetTimes(0, t_end);
        solver.SetTimeStep(0.001);

        solver.SetInitialCondition(initial_condition);

        Vec result = solver.Solve();
        ReplicatableVector result_repl(result);

        // Check solution is u = e^{-3*t*pi*pi} sin(x*pi)*sin(y*pi)*sin(z*pi), t=0.1
        for (unsigned i=0; i<result_repl.GetSize(); i++)
        {
            double x = mesh.GetNode(i)->GetPoint()[0];
            double y = mesh.GetNode(i)->GetPoint()[1];
            double z = mesh.GetNode(i)->GetPoint()[2];
            double u = exp(-3*t_end*M_PI*M_PI)*sin(x*M_PI)*sin(y*M_PI)*sin(z*M_PI);
            TS_ASSERT_DELTA(result_repl[i], u, 0.1);
        }

        PetscTools::Destroy(initial_condition);
        PetscTools::Destroy(result);
    }

    void TestProblemChecksUsingBathWithMultipleBathConductivities()
    {
        TrianglesMeshReader<2,2> reader("mesh/test/data/2D_0_to_1mm_400_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(reader);

        std::set<unsigned> tissue_ids;
        tissue_ids.insert(0);

        std::set<unsigned> bath_ids;
        bath_ids.insert(1);
        bath_ids.insert(2); // non-default identifier!

        HeartConfig::Instance()->SetTissueAndBathIdentifiers(tissue_ids, bath_ids);

        BathCellFactory<2> cell_factory( 0.0, Create_c_vector(0.0, 0.0) );
        BidomainProblem<2> bidomain_problem( &cell_factory ); // non-bath problem, despite specifying bath stuff above!
        bidomain_problem.SetMesh( &mesh );
        TS_ASSERT_THROWS_THIS( bidomain_problem.Initialise() , "User has set bath identifiers, but the BidomainProblem isn't expecting a bath. Did you mean to use BidomainProblem(..., true)? Or alternatively, BidomainWithBathProblem(...)?");
    }

std::vector<unsigned> NonlinearElasticityTools<DIM>::GetNodesByComponentValue(TetrahedralMesh<DIM,DIM>& rMesh,
                                                          unsigned component,
                                                          double value)
{
    std::vector<unsigned> fixed_nodes;
    double tol = 1e-8;
    for (unsigned i=0; i<rMesh.GetNumNodes(); i++)
    {
        if ( fabs(rMesh.GetNode(i)->rGetLocation()[component] - value)<1e-8)
        {
            fixed_nodes.push_back(i);
        }
    }

    if (fixed_nodes.size() == 0)
    {
        EXCEPTION("Could not find any nodes on requested surface (note: tolerance = "<<tol<<")");
    }

    return fixed_nodes;
}

    void TestDefineZeroDirichletOnMeshBoundary()
    {
        // Load a 2D square mesh with 1 central non-boundary node
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_4_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        BoundaryConditionsContainer<2,2,1> bcc;

        bcc.DefineZeroDirichletOnMeshBoundary(&mesh);

        // Check boundary nodes have the right condition
        for (int i=0; i<4; i++)
        {
            double value = bcc.GetDirichletBCValue(mesh.GetNode(i));
            TS_ASSERT_DELTA(value, 0.0, 1e-12);
        }

        // Check non-boundary node has no condition
        TS_ASSERT(!bcc.HasDirichletBoundaryCondition(mesh.GetNode(4)));
    }

    // This test covers the case when the hdf5 file contains more than 3 variables
    void TestMeshalyzerConversionLotsOfVariables() throw(Exception)
    {
        std::string output_dir = "TestHdf5Converters_TestMeshalyzerConversionLotsOfVariables";

        /*
         * Firstly, copy the .h5 file to CHASTE_TEST_OUTPUT/TestHdf5ToMeshalyzerConverter,
         * as that is where the reader reads from.
         */
        CopyToTestOutputDirectory("heart/test/data/many_variables/many_variables.h5", output_dir);

        TrianglesMeshReader<1,1> mesh_reader("heart/test/data/many_variables/1D_65_elements");
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        // Convert
        Hdf5ToMeshalyzerConverter<1,1> converter(FileFinder(output_dir, RelativeTo::ChasteTestOutput),
                                                 "many_variables", &mesh, true);

        std::vector<std::string> variable_names;
        variable_names.push_back("V");
        variable_names.push_back("I_ks");
        variable_names.push_back("I_kr");
        variable_names.push_back("I_Ca_tot");
        variable_names.push_back("I_tot");
        variable_names.push_back("I_Na_tot");

        std::string test_output_directory = OutputFileHandler::GetChasteTestOutputDirectory();
        for (unsigned i=0; i<variable_names.size(); i++)
        {
            // Compare the results files
            FileComparison(test_output_directory + "/" + output_dir + "/output/many_variables_"
                           + variable_names[i] + ".dat",
                           "heart/test/data/many_variables/many_variables_"
                           + variable_names[i] + ".dat").CompareFiles();
        }

        // Compare the time information file
        FileComparison(test_output_directory + output_dir + "/output/many_variables_times.info",
                       "heart/test/data/many_variables/many_variables_times.info").CompareFiles();
    }

    void Test2DMeshRotation()
    {
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/2D_0_to_1mm_200_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        double angle = M_PI;

        mesh.Rotate(angle);

        TetrahedralMesh<2,2> original_mesh;
        original_mesh.ConstructFromMeshReader(mesh_reader);

        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            // Find new coordinates of the translated node
            Node<2>* p_node = mesh.GetNode(i);
            ChastePoint<2> new_coordinate = p_node->GetPoint();

            // Get original node
            Node<2>* p_original_node = original_mesh.GetNode(i);
            ChastePoint<2> original_coordinate = p_original_node->GetPoint();

            // Run a test to make sure the node has gone to the correct place
            TS_ASSERT_DELTA(original_coordinate[0], -new_coordinate[0], 1e-5);
            TS_ASSERT_DELTA(original_coordinate[1], -new_coordinate[1], 1e-5);
        }

        // Check volume conservation
        double mesh_volume = mesh.GetVolume();
        double original_mesh_volume = original_mesh.GetVolume();

        TS_ASSERT_DELTA(mesh_volume, original_mesh_volume, 1e-5);
    }

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

        // Pick a random node and store spatial position
        Node<3>* p_node = mesh.GetNode(10);
        ChastePoint<3> original_coordinate = p_node->GetPoint();

        double mesh_volume = mesh.GetVolume();

        const double x_movement = 1.0;
        const double y_movement = 2.5;
        const double z_movement = -3.75;
        mesh.Translate(x_movement, y_movement, z_movement);

        ChastePoint<3>  new_coordinate = p_node->GetPoint();
        double new_mesh_volume = mesh.GetVolume();

        TS_ASSERT_DELTA(mesh_volume, new_mesh_volume, 1e-6);
        TS_ASSERT_DELTA(original_coordinate[0], new_coordinate[0]-x_movement, 1e-6);
        TS_ASSERT_DELTA(original_coordinate[1], new_coordinate[1]-y_movement, 1e-6);
        TS_ASSERT_DELTA(original_coordinate[2], new_coordinate[2]-z_movement, 1e-6);
    }

    void TestSimpleOrthotropicNotDistributed() throw (Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/box_shaped_heart/box_heart");
        std::string epi_face_file = "heart/test/data/box_shaped_heart/epi.tri";
        std::string rv_face_file = "heart/test/data/box_shaped_heart/rv.tri";
        std::string lv_face_file = "heart/test/data/box_shaped_heart/lv.tri";

        TetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        StreeterFibreGenerator<3> fibre_generator(mesh);
        fibre_generator.SetSurfaceFiles(epi_face_file, rv_face_file, lv_face_file, false);
        fibre_generator.SetApexToBase(0);

        OutputFileHandler handler("shorter_streeter", false);
        fibre_generator.WriteData(handler, "box_heart_not_dist.ortho");

        FileFinder fibre_file1 = handler.FindFile("box_heart_not_dist.ortho");
        FileFinder fibre_file2("heart/test/data/box_shaped_heart/box_heart.ortho", RelativeTo::ChasteSourceRoot);

        CompareGeneratedWithReferenceFile(fibre_file1, ORTHO, fibre_file2, ORTHO);
    }

    void TestMonodomainMeshalyzerConversion() throw(Exception)
    {
        // Firstly, copy ./heart/test/data/MonoDg01d/*.h5 to CHASTE_TEST_OUTPUT/TestHdf5ToMeshalyzerConverter,
        // as that is where the reader reads from.
        std::string output_folder("TestHdf5Converters_TestMonodomainMeshalyzerConversion");
        CopyToTestOutputDirectory("heart/test/data/Monodomain1d/MonodomainLR91_1d.h5", output_folder);

        TrianglesMeshReader<1,1> mesh_reader("mesh/test/data/1D_0_to_1_100_elements");
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        // Convert
        Hdf5ToMeshalyzerConverter<1,1> converter(FileFinder(output_folder,RelativeTo::ChasteTestOutput),
                                                 "MonodomainLR91_1d", &mesh, true, 10 /* precision specified for coverage */);

        // Compare the voltage file with a correct version
        std::string test_output_directory = OutputFileHandler::GetChasteTestOutputDirectory();
        FileComparison(test_output_directory + output_folder + "/output/MonodomainLR91_1d_V.dat",
                       "heart/test/data/Monodomain1d/MonodomainLR91_1d_V.dat").CompareFiles();
        FileComparison(test_output_directory + output_folder + "/output/MonodomainLR91_1d_times.info",
                       "heart/test/data/Monodomain1d/MonodomainLR91_1d_times.info").CompareFiles();
    }

    void TestRemeshFullTree() throw(Exception)
    {
        TrianglesMeshReader<1,3> reader("lung/test/data/TestSubject002");

        TetrahedralMesh<1,3> mesh;
        mesh.ConstructFromMeshReader(reader);

        AirwayPropertiesCalculator calculator(mesh, 0);

        TS_ASSERT_EQUALS( mesh.GetNumNodes(), 136625u);
        TS_ASSERT_EQUALS( mesh.GetNumElements(), 136624u);

        //Create remesher object
        AirwayRemesher remesher(mesh, 0u);

        MutableMesh<1,3> output_mesh_one;
        remesher.Remesh(output_mesh_one, calculator.GetBranches()[0]->GetPoiseuilleResistance()*1e7); //Key the tolerance relative to the trachea

        TS_ASSERT_EQUALS( output_mesh_one.GetNumNodes(), 168045u);
        TS_ASSERT_EQUALS( output_mesh_one.GetNumElements(), 168044u);

//        //To visualise
//
//        VtkMeshWriter<1,3> writer("TestAirwayRemesher", "Novartis002_remeshed");
//        std::vector<double> radii(output_mesh_one.GetNumElements());
//
//        for(TetrahedralMesh<1,3>::ElementIterator iter = output_mesh_one.GetElementIteratorBegin();
//            iter != output_mesh_one.GetElementIteratorEnd();
//            ++iter)
//        {
//            radii[iter->GetIndex()] = iter->GetAttribute();
//        }
//
//        writer.AddCellData("radii", radii);
//        writer.WriteFilesUsingMesh(output_mesh_one);
//
//        TrianglesMeshWriter<1,3> writer2("TestAirwayRemesher", "Novartis002_remeshed", false);
//        writer2.WriteFilesUsingMesh(output_mesh_one);
    }