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

    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 TestCellwiseDataGradientVerySmallMesh() throw(Exception)
    {
        // Create a simple mesh
        TrianglesMeshReader<2,2> mesh_reader("mesh/test/data/square_2_elements");
        MutableMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

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

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

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

        // With the algorithm being used, the numerical gradient is (1,0)
        // for each of the nodes
        for (unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            TS_ASSERT_DELTA(gradient.rGetGradient(i)(0), 1.0, 1e-9);
            TS_ASSERT_DELTA(gradient.rGetGradient(i)(1), 0.0, 1e-9);
        }
    }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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