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

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

        double mesh_volume = mesh.GetVolume();

        mesh.Scale(1.0);
        TS_ASSERT_DELTA(mesh_volume, mesh.GetVolume(), 1e-6);

        mesh.Scale(2.0, 3.0, 4.0);
        TS_ASSERT_DELTA(24.0*mesh_volume, mesh.GetVolume(), 1e-6);

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

        mesh.Scale(0.5, 1.0/3.0, 0.25);

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

        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 TestValidate()
    {
        // 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;

        // No BCs yet, so shouldn't validate
        TS_ASSERT(!bcc.Validate(&mesh));

        // Add some BCs
        ConstBoundaryCondition<2> *bc = new ConstBoundaryCondition<2>(0.0);
        bcc.AddDirichletBoundaryCondition(mesh.GetNode(0), bc);
        bcc.AddDirichletBoundaryCondition(mesh.GetNode(1), bc);
        bcc.AddDirichletBoundaryCondition(mesh.GetNode(3), bc);
        TetrahedralMesh<2,2>::BoundaryElementIterator iter
        = mesh.GetBoundaryElementIteratorEnd();
        iter--;
        bcc.AddNeumannBoundaryCondition(*iter, bc); // 2 to 3
        iter--;
        bcc.AddNeumannBoundaryCondition(*iter, bc); // 1 to 2

        TS_ASSERT(bcc.Validate(&mesh));
    }

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

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

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

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

        // Translate
        mesh.Translate(trans_vec);
        c_vector<double, 3> new_location1 = point1.rGetLocation();
        c_vector<double, 3> 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 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);
    }

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

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