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

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

    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 TestBidomainProblemWithDistributedMeshFromMemfem3DParMetis() throw(Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(1);  //ms
        HeartConfig::Instance()->SetOutputDirectory("DistributedMesh3dRepViaTri");
        HeartConfig::Instance()->SetOutputFilenamePrefix("tetrahedral3d");

        // The default stimulus in PlaneStimulusCellFactory is not enough to generate propagation
        // here, increasing it an order of magnitude
        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 3> cell_factory(-6000);

        // To avoid an issue with the Event handler only one simulation should be
        // in existance at a time: therefore monodomain simulation is defined in a block
        double seq_ave_voltage=0.0;
        {
            ///////////////////////////////////////////////////////////////////
            // TetrahedralMesh from Triangles
            ///////////////////////////////////////////////////////////////////
            TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/SlabFromMemfem");
            TetrahedralMesh<3,3> mesh;

            mesh.ConstructFromMeshReader(mesh_reader);
            TS_ASSERT_DELTA(mesh.GetVolume(), 1.5625, 1e-6);
            TS_ASSERT_EQUALS(mesh_reader.GetNumNodes(), 381u);
            TS_ASSERT_EQUALS(mesh_reader.GetNumElements(), 1030u);
            TS_ASSERT_EQUALS(mesh_reader.GetNumFaces(), 758u);


            BidomainProblem<3> nondistributed_problem( &cell_factory );
            nondistributed_problem.SetMesh(&mesh);
            nondistributed_problem.Initialise();

            HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
            HeartConfig::Instance()->SetCapacitance(1.0);

            nondistributed_problem.Solve();

            DistributedVector dist_nondistributed_voltage = nondistributed_problem.GetSolutionDistributedVector();
            DistributedVector::Stripe nondistributed_voltage(dist_nondistributed_voltage, 0);
            DistributedVector::Stripe nondistributed_potential(dist_nondistributed_voltage, 1);

            double seq_local_ave_voltage = 0.0;

            for (DistributedVector::Iterator index = dist_nondistributed_voltage.Begin();
                 index != dist_nondistributed_voltage.End();
                 ++index)
            {
                if (index.Global==0)
                {
                    TS_ASSERT_LESS_THAN(0, nondistributed_voltage[index]);
                }

                seq_local_ave_voltage += nondistributed_voltage[index];
            }

            MPI_Reduce(&seq_local_ave_voltage, &seq_ave_voltage, 1, MPI_DOUBLE, MPI_SUM, PetscTools::MASTER_RANK, PETSC_COMM_WORLD);
            seq_ave_voltage /= mesh.GetNumNodes();
        }

        ///////////////////////////////////////////////////////////////////
        // DistributedTetrahedralMesh from Memfem
        ///////////////////////////////////////////////////////////////////
        HeartConfig::Instance()->SetOutputDirectory("DistributedMesh3dDistViaMem");

        BidomainProblem<3> distributed_problem( &cell_factory );

        HeartConfig::Instance()->SetMeshFileName("mesh/test/data/Memfem_slab");
        distributed_problem.Initialise();

        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        HeartConfig::Instance()->SetCapacitance(1.0);

        distributed_problem.Solve();

        DistributedVector dist_distributed_voltage = distributed_problem.GetSolutionDistributedVector();
        DistributedVector::Stripe distributed_voltage(dist_distributed_voltage, 0);
        DistributedVector::Stripe distributed_potential(dist_distributed_voltage, 1);

        double para_local_ave_voltage = 0.0;

        for (DistributedVector::Iterator index = dist_distributed_voltage.Begin();
             index != dist_distributed_voltage.End();
             ++index)
        {
            if (index.Global==0)
            {
                TS_ASSERT_LESS_THAN(0, distributed_voltage[index]);
            }

            para_local_ave_voltage += distributed_voltage[index];
        }


        double para_ave_voltage;
        MPI_Reduce(&para_local_ave_voltage, &para_ave_voltage, 1, MPI_DOUBLE, MPI_SUM, PetscTools::MASTER_RANK, PETSC_COMM_WORLD);
        para_ave_voltage /= distributed_problem.rGetMesh().GetNumNodes();

        ///////////////////////////////////////////////////////////////////
        // compare
        ///////////////////////////////////////////////////////////////////
        if (PetscTools::AmMaster())
//.........这里部分代码省略.........