C++ TetrahedralMesh::GetNumNodes方法代码示例

    void TestGenerateBasicRandomWithNoSpecifiedProliferativeCellType() throw(Exception)
    {
        // Create mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_2_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

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

        // Test that cells were generated correctly
        TS_ASSERT_EQUALS(cells.size(), mesh.GetNumNodes());

        for (unsigned i=0; i<cells.size(); i++)
        {
            TS_ASSERT_EQUALS(cells[i]->GetCellCycleModel()->GetDimension(), 2u);
            TS_ASSERT_EQUALS(cells[i]->GetCellProliferativeType()->IsType<StemCellProliferativeType>(), true);
            // Should lie between -24 and 0
            double birth_time=cells[i]->GetBirthTime();
            ///\todo Breaks Intel 10? TS_ASSERT_LESS_THAN_EQUALS(birth_time, 0.0);
            TS_ASSERT_LESS_THAN_EQUALS(-24.0, birth_time);
        }
    }

    void TestGenerateBasicRandomWithFixedDurationGenerationBasedCellCycleModel() throw(Exception)
    {
        // Create mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_2_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

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

        // Test that cells were generated correctly
        TS_ASSERT_EQUALS(cells.size(), mesh.GetNumNodes());

        for (unsigned i=0; i<cells.size(); i++)
        {
            // Should lie between -24 and 0
            TS_ASSERT_LESS_THAN_EQUALS(cells[i]->GetBirthTime(), 0.0);
            TS_ASSERT_LESS_THAN_EQUALS(-24.0, cells[i]->GetBirthTime());
            TS_ASSERT_EQUALS(cells[i]->GetCellCycleModel()->GetDimension(), 2u);
            TS_ASSERT_EQUALS(cells[i]->GetCellProliferativeType(), p_transit_type);
        }

        // Test exact random numbers as test re-seeds random number generator.
        TS_ASSERT_DELTA(cells[0]->GetBirthTime(), -7.1141, 1e-4);
        TS_ASSERT_DELTA(cells[1]->GetBirthTime(), -10.1311, 1e-4);
        TS_ASSERT_DELTA(cells[2]->GetBirthTime(), -10.2953, 1e-4);
    }

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

    /**
     * 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 TestDistancesToFaceDumb()
    {
        TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_21_nodes_side/Cube21"); // 5x5x5mm cube (internode distance = 0.25mm)

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

        TS_ASSERT_EQUALS(mesh.GetNumNodes(), 9261u); // 21x21x21 nodes
        TS_ASSERT_EQUALS(mesh.GetNumElements(), 48000u);
        TS_ASSERT_EQUALS(mesh.GetNumBoundaryElements(), 4800u);

        DistributedTetrahedralMesh<3,3> parallel_mesh(DistributedTetrahedralMeshPartitionType::DUMB); // No reordering
        parallel_mesh.ConstructFromMeshReader(mesh_reader);
        TS_ASSERT_EQUALS(parallel_mesh.GetNumNodes(), 9261u); // 21x21x21 nodes
        TS_ASSERT_EQUALS(parallel_mesh.GetNumElements(), 48000u);
        TS_ASSERT_EQUALS(parallel_mesh.GetNumBoundaryElements(), 4800u);


        std::vector<unsigned> map_left;
        for (unsigned index=0; index<mesh.GetNumNodes(); index++)
        {
            // Get the nodes at the left face of the cube
            if (mesh.GetNode(index)->rGetLocation()[0] + 0.25 < 1e-6)
            {
                map_left.push_back(index);
            }
        }

        TS_ASSERT_EQUALS(map_left.size(), 21u*21u);

        DistanceMapCalculator<3,3> distance_calculator(mesh);
        std::vector<double> distances;
        distance_calculator.ComputeDistanceMap(map_left, distances);

        DistanceMapCalculator<3,3> parallel_distance_calculator(parallel_mesh);
        std::vector<double> parallel_distances;
        parallel_distance_calculator.ComputeDistanceMap(map_left, parallel_distances);

        TS_ASSERT_EQUALS(distance_calculator.mRoundCounter, 1u);
        TS_ASSERT_DELTA(parallel_distance_calculator.mRoundCounter, 2u, 1u);// 1 2 or 3

        for (unsigned index=0; index<distances.size(); index++)
        {
            // The distance should be equal to the x-coordinate of the point (minus the offset of the left face of the cube)
            c_vector<double, 3> node = mesh.GetNode(index)->rGetLocation();
            TS_ASSERT_DELTA(distances[index], node[0]+0.25,1e-11);
            TS_ASSERT_DELTA(parallel_distances[index], node[0]+0.25,1e-11);
        }
    }

    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 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 TestXaxisRotation3DWithHomogeneousUblas()
    {
        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

        c_matrix<double, 4, 4> x_rotation_matrix = identity_matrix<double>(4);

        double theta = M_PI/2;

        x_rotation_matrix(1,1) = cos(theta);
        x_rotation_matrix(1,2) = sin(theta);
        x_rotation_matrix(2,1) = -sin(theta);
        x_rotation_matrix(2,2) = cos(theta);

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

        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            Node<3>* p_node = mesh.GetNode(i);
            ChastePoint<3> point = p_node->GetPoint();

            c_vector<double, 4> point_location;

            point_location[0] = point[0];
            point_location[1] = point[1];
            point_location[2] = point[2];
            point_location[3] = 1.0;

            c_vector<double, 4> new_point_location = prod(x_rotation_matrix, point_location);

            TS_ASSERT_EQUALS(new_point_location[3], 1.0);

            point.SetCoordinate(0, new_point_location[0]);
            point.SetCoordinate(1, new_point_location[1]);
            point.SetCoordinate(2, new_point_location[2]);
            p_node->SetPoint(point);
        }

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

        mesh.RefreshMesh();

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

    void TestReadMeshes(void) throw(Exception)
    {
        {
            READER_2D reader("mesh/test/data/square_4_elements_gmsh.msh");
            TetrahedralMesh<2,2> mesh;
            mesh.ConstructFromMeshReader(reader);
            TS_ASSERT_EQUALS(mesh.GetNumNodes(), 5u);
            TS_ASSERT_EQUALS(mesh.GetNumElements(), 4u);
            TS_ASSERT_EQUALS(mesh.GetNumBoundaryElements(), 4u);
        }

        {
            READER_3D reader("mesh/test/data/simple_cube_gmsh.msh");
            TetrahedralMesh<3,3> mesh;
            mesh.ConstructFromMeshReader(reader);
            TS_ASSERT_EQUALS(mesh.GetNumNodes(), 14u);
            TS_ASSERT_EQUALS(mesh.GetNumElements(), 24u);
            TS_ASSERT_EQUALS(mesh.GetNumBoundaryElements(), 24u);
        }

        {
           READER_2D reader("mesh/test/data/quad_square_4_elements_gmsh.msh",2,2);
           QuadraticMesh<2> mesh;
           mesh.ConstructFromMeshReader(reader);
           TS_ASSERT_EQUALS(mesh.GetNumNodes(), 13u);
           TS_ASSERT_EQUALS(mesh.GetNumElements(), 4u);
           TS_ASSERT_EQUALS(mesh.GetNumBoundaryElements(), 4u);
        }

        {
           READER_3D reader("mesh/test/data/quad_cube_gmsh.msh",2,2);
           QuadraticMesh<3> mesh;
           mesh.ConstructFromMeshReader(reader);
           TS_ASSERT_EQUALS(mesh.GetNumNodes(), 63u);
           TS_ASSERT_EQUALS(mesh.GetNumElements(), 24u);
           TS_ASSERT_EQUALS(mesh.GetNumBoundaryElements(), 24u);
        }
    }

    void Test3dBathIntracellularStimulation()
    {
        HeartConfig::Instance()->SetSimulationDuration(1);  //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainBath3d");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain_bath_3d");

        c_vector<double,3> centre;
        centre(0) = 0.05;
        centre(1) = 0.05;
        centre(2) = 0.05;
        BathCellFactory<3> cell_factory(-2.5e7, centre); // stimulates x=0.05 node

        BidomainProblem<3> bidomain_problem( &cell_factory, true );

        TetrahedralMesh<3,3> mesh;
        mesh.ConstructRegularSlabMesh(0.01, 0.1, 0.1, 0.1);

        // Set everything outside a central sphere (radius 0.4) to be bath
        for (unsigned i=0; i<mesh.GetNumElements(); i++)
        {
            double x = mesh.GetElement(i)->CalculateCentroid()[0];
            double y = mesh.GetElement(i)->CalculateCentroid()[1];
            double z = mesh.GetElement(i)->CalculateCentroid()[2];
            if (sqrt((x-0.05)*(x-0.05) + (y-0.05)*(y-0.05) + (z-0.05)*(z-0.05)) > 0.04)
            {
                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);
            }
        }

        // a hardcoded value
        TS_ASSERT_DELTA(sol_repl[2*404], 39.6833, 1e-3);
    }

    void TestGenerateBasicWithFixedDurationGenerationBasedCellCycleModel() throw(Exception)
    {
        // Create mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_2_elements");
        TetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

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

        // Test that cells were generated correctly
        TS_ASSERT_EQUALS(cells.size(), mesh.GetNumNodes());

        for (unsigned i=0; i<cells.size(); i++)
        {
            TS_ASSERT_DELTA(cells[i]->GetBirthTime(), -(double)(i), 1e-9);
            TS_ASSERT_EQUALS(cells[i]->GetCellCycleModel()->GetDimension(), 2u);
        }

        // Test with extra input argument
        std::vector<unsigned> location_indices;
        location_indices.push_back(2);
        location_indices.push_back(7);
        location_indices.push_back(9);

        std::vector<CellPtr> cells2;
        CellsGenerator<FixedDurationGenerationBasedCellCycleModel, 2> cells_generator2;
        cells_generator2.GenerateBasic(cells2, 3, location_indices);

        TS_ASSERT_EQUALS(cells2.size(), 3u);
        TS_ASSERT_DELTA(cells2[0]->GetBirthTime(), -2.0, 1e-4);
        TS_ASSERT_DELTA(cells2[1]->GetBirthTime(), -7.0, 1e-4);
        TS_ASSERT_DELTA(cells2[2]->GetBirthTime(), -9.0, 1e-4);
    }

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

        // Translations - add a constant vector to each node

        c_vector<double,3>  displacement;
        displacement(0) = 1;
        displacement(1) = 1;
        displacement(2) = 1;

        // Translate the mesh along the vector displacement
        mesh.Translate(displacement);

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

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

            // Get original node
            Node<3>* p_original_node = original_mesh.GetNode(i);
            ChastePoint<3> 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] + displacement[0], new_coordinate[0], 1e-5);
            TS_ASSERT_DELTA(original_coordinate[1] + displacement[1], new_coordinate[1], 1e-5);
            TS_ASSERT_DELTA(original_coordinate[2] + displacement[2], new_coordinate[2], 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);
    }

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 Test3DAngleAxisRotation()
    {
        TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_136_elements");
        TetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        c_vector<double,3> axis;
        axis(0) = 1;
        axis(1) = 0;
        axis(2) = 0;

        double angle = M_PI;

        mesh.Rotate(axis, angle);

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

        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            Node<3>* p_node = mesh.GetNode(i);
            ChastePoint<3> new_coordinate = p_node->GetPoint();

            // Get original node
            Node<3>* p_original_node = original_mesh.GetNode(i);
            ChastePoint<3>  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);
            TS_ASSERT_DELTA(original_coordinate[2], -new_coordinate[2], 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 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);
    }