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

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

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

    void TestDistancesToCorner() throw (Exception)
    {
        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);

        unsigned num_nodes=9261u;
        TS_ASSERT_EQUALS(mesh.GetNumNodes(), num_nodes); // 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(), num_nodes); // 21x21x21 nodes
        TS_ASSERT_EQUALS(parallel_mesh.GetNumElements(), 48000u);
        TS_ASSERT_EQUALS(parallel_mesh.GetNumBoundaryElements(), 4800u);

        unsigned far_index=9260u;
        c_vector<double,3> far_corner=mesh.GetNode(far_index)->rGetLocation();
        TS_ASSERT_DELTA( far_corner[0], 0.25, 1e-11);
        TS_ASSERT_DELTA( far_corner[1], 0.25, 1e-11);
        TS_ASSERT_DELTA( far_corner[2], 0.25, 1e-11);
        try
        {
            c_vector<double,3> parallel_far_corner=parallel_mesh.GetNode(far_index)->rGetLocation();
            TS_ASSERT_DELTA( parallel_far_corner[0], 0.25, 1e-11);
            TS_ASSERT_DELTA( parallel_far_corner[1], 0.25, 1e-11);
            TS_ASSERT_DELTA( parallel_far_corner[2], 0.25, 1e-11);
        }
        catch (Exception&)
        {
        }

        std::vector<unsigned> map_far_corner;
        map_far_corner.push_back(far_index);

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

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

        TS_ASSERT_EQUALS(distance_calculator.mRoundCounter, 1u);
        //Nodes in mesh are order such that a dumb partitioning will give a sequential handover from proc0 to proc1...
        TS_ASSERT_EQUALS(parallel_distance_calculator.mRoundCounter, PetscTools::GetNumProcs());
        //Note unsigned division is okay here
        TS_ASSERT_DELTA(parallel_distance_calculator.mPopCounter, num_nodes/PetscTools::GetNumProcs(), 1u);
        TS_ASSERT_DELTA(distance_calculator.mPopCounter, num_nodes, 1u);
        for (unsigned index=0; index<distances.size(); index++)
        {
            c_vector<double, 3> node = mesh.GetNode(index)->rGetLocation();

            //Straightline distance
            double euclidean_distance = norm_2(far_corner - node);
            // x + y + z distance
            double manhattan_distance = norm_1(far_corner - node);
            //If they differ, then allow the in-mesh distance to be in between
            double error_bound = (manhattan_distance - euclidean_distance)/2.0;
            //If they don't differ, then we expect the in-mesh distance to be similar
            if (error_bound < 1e-15)
            {
                error_bound = 1e-15;
            }
            TS_ASSERT_LESS_THAN_EQUALS(distances[index], manhattan_distance+DBL_EPSILON);
            TS_ASSERT_LESS_THAN_EQUALS(euclidean_distance, distances[index]+DBL_EPSILON);
            TS_ASSERT_DELTA(distances[index], euclidean_distance, error_bound);

            TS_ASSERT_DELTA(distances[index], parallel_distances[index], 1e-15);
        }

        // Test some point-to-point distances
        RandomNumberGenerator::Instance()->Reseed(1);
        unsigned trials=25;
        unsigned pops=0;
        unsigned sequential_pops=0;
        for (unsigned i=0; i<trials; i++)
        {
            unsigned index=RandomNumberGenerator::Instance()->randMod(parallel_distances.size());
            TS_ASSERT_DELTA(parallel_distance_calculator.SingleDistance(9260u, index), parallel_distances[index], 1e-15);
            TS_ASSERT_DELTA(distance_calculator.SingleDistance(9260u, index), parallel_distances[index], 1e-15);
            pops += parallel_distance_calculator.mPopCounter;
            sequential_pops += distance_calculator.mPopCounter;
            TS_ASSERT_LESS_THAN_EQUALS(parallel_distance_calculator.mRoundCounter, PetscTools::GetNumProcs()+2);
        }

        // Without A*: TS_ASSERT_DELTA(sequential_pops/(double)trials, num_nodes/2, 300);
        TS_ASSERT_LESS_THAN(sequential_pops/(double)trials, num_nodes/20.0);
        if (PetscTools::IsSequential())
        {
            //Early termination
            TS_ASSERT_EQUALS(pops, sequential_pops);
        }
        else
        {
            //Early termination on remote processes is not yet possible
            //This may lead to multiple updates from remote
            //A* Leads to even more updates on average
//.........这里部分代码省略.........

    void TestGetRadiusVectorsAndConstructStructureTensors(void) throw(Exception)
    {
        // Test in three parts to use the results of one test in the next...
        //
        //

        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

        calculator.GetRadiusVectors();

        std::vector< c_vector<double,3> >& radius_vectors = calculator.mRadiusVectors;

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

        TS_ASSERT_DELTA(radius_vectors[5][0], 0.475, 1e-9);
        TS_ASSERT_DELTA(radius_vectors[5][1], 0.225, 1e-9);
        TS_ASSERT_DELTA(radius_vectors[5][2], 0.475, 1e-9);

        TS_ASSERT_EQUALS(radius_vectors.size(), mesh.GetNumElements());

        ///////////////////////////////////////////////////////////
        // Test ConstructStructureTensors()
        ///////////////////////////////////////////////////////////
        calculator.ConstructStructureTensors();
        std::vector< c_matrix<double,3,3> >& tensor_i = calculator.mStructureTensors;

        // Worked out by hand...
        TS_ASSERT_DELTA(tensor_i[0](0,0),7.5625e-02,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](0,1),6.8750e-03,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](0,2),7.5625e-02,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](1,0),6.8750e-03,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](1,1),6.2500e-04,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](1,2),6.8750e-03,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](2,0),7.5625e-02,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](2,1),6.8750e-03,1e-9);
        TS_ASSERT_DELTA(tensor_i[0](2,2),7.5625e-02,1e-9);

        TS_ASSERT_DELTA(tensor_i[5](0,0),0.225625,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](0,1),0.106875,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](0,2),0.225625,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](1,0),0.106875,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](1,1),0.050625,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](1,2),0.106875,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](2,0),0.225625,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](2,1),0.106875,1e-9);
        TS_ASSERT_DELTA(tensor_i[5](2,2),0.225625,1e-9);

        //////////////////////////////////////////////////////////////
        // Test SmoothStructureTensor()
        //////////////////////////////////////////////////////////////
        calculator.SmoothStructureTensors();
        std::vector< c_matrix<double,3,3> >& tensor_smooth = calculator.mSmoothedStructureTensors;

        // hard-coded (as difficult to test)
        TS_ASSERT_DELTA(tensor_smooth[0](0,0), 0.075625, 1e-5);
    }

    // In this test we have no cardiac tissue, so that the equations are just sigma * phi_e''=0
    // throughout the domain (with a Neumann boundary condition on x=1 and a dirichlet boundary
    // condition (ie grounding) on x=0), so the exact solution can be calculated and compared
    // against.
    void Test1dProblemOnlyBathGroundedOneSide() throw (Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(0.5);  //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainBathOnlyBath");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain_bath");

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

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

        for(unsigned i=0; i<mesh.GetNumElements(); i++)
        {
            mesh.GetElement(i)->SetAttribute(HeartRegionCode::GetValidBathId());
        }

        // create boundary conditions container
        double boundary_val = 1.0;
        boost::shared_ptr<BoundaryConditionsContainer<1,1,2> > p_bcc(new BoundaryConditionsContainer<1,1,2>);
        ConstBoundaryCondition<1>* p_bc_stim = new ConstBoundaryCondition<1>(boundary_val);
        ConstBoundaryCondition<1>* p_zero_stim = new ConstBoundaryCondition<1>(0.0);

        // loop over boundary elements and set (sigma\gradphi).n = 1.0 on RHS edge
        for(TetrahedralMesh<1,1>::BoundaryElementIterator iter
              = mesh.GetBoundaryElementIteratorBegin();
           iter != mesh.GetBoundaryElementIteratorEnd();
           iter++)
        {
            if (((*iter)->GetNodeLocation(0))[0]==1.0)
            {
                /// \todo: I think you need to provide a boundary condition for unknown#1 if you are gonig to provide one for unknown#2?
                p_bcc->AddNeumannBoundaryCondition(*iter, p_zero_stim, 0);
                p_bcc->AddNeumannBoundaryCondition(*iter, p_bc_stim,   1);
            }
        }

        BidomainWithBathProblem<1> bidomain_problem( &cell_factory );

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

        // fix phi=0 on LHS edge
        std::vector<unsigned> fixed_nodes;
        fixed_nodes.push_back(0);
        bidomain_problem.SetFixedExtracellularPotentialNodes(fixed_nodes);

        bidomain_problem.Solve();

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

        // test phi = x*boundary_val/sigma (solution of phi''=0, phi(0)=0, sigma*phi'(1)=boundary_val
        for(unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            double bath_cond = HeartConfig::Instance()->GetBathConductivity();
            double x = mesh.GetNode(i)->rGetLocation()[0];
            TS_ASSERT_DELTA(sol_repl[2*i],   0.0,   1e-12);               // V
            TS_ASSERT_DELTA(sol_repl[2*i+1], x*boundary_val/bath_cond, 1e-4);   // phi_e
        }
    }