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

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

VoltageInterpolaterOntoMechanicsMesh<DIM>::VoltageInterpolaterOntoMechanicsMesh(
                                     TetrahedralMesh<DIM,DIM>& rElectricsMesh,
                                     QuadraticMesh<DIM>& rMechanicsMesh,
                                     std::vector<std::string>& rVariableNames,
                                     std::string directory,
                                     std::string inputFileNamePrefix)
{
    // Read the data from the HDF5 file
    Hdf5DataReader reader(directory,inputFileNamePrefix);

    unsigned num_timesteps = reader.GetUnlimitedDimensionValues().size();

    // set up the elements and weights for the coarse nodes in the fine mesh
    FineCoarseMeshPair<DIM> mesh_pair(rElectricsMesh, rMechanicsMesh);
    mesh_pair.SetUpBoxesOnFineMesh();
    mesh_pair.ComputeFineElementsAndWeightsForCoarseNodes(true);
    assert(mesh_pair.rGetElementsAndWeights().size()==rMechanicsMesh.GetNumNodes());

    // create and setup a writer
    Hdf5DataWriter* p_writer = new Hdf5DataWriter(*rMechanicsMesh.GetDistributedVectorFactory(),
                                                  directory,
                                                  "voltage_mechanics_mesh",
                                                  false, //don't clean
                                                  false);

    std::vector<int> columns_id;
    for (unsigned var_index = 0; var_index < rVariableNames.size(); var_index++)
    {
        std::string var_name = rVariableNames[var_index];
        columns_id.push_back( p_writer->DefineVariable(var_name,"mV") );
    }

    p_writer->DefineUnlimitedDimension("Time","msecs", num_timesteps);
    p_writer->DefineFixedDimension( rMechanicsMesh.GetNumNodes() );
    p_writer->EndDefineMode();

    assert(columns_id.size() == rVariableNames.size());

    // set up a vector to read into
    DistributedVectorFactory factory(rElectricsMesh.GetNumNodes());
    Vec voltage = factory.CreateVec();
    std::vector<double> interpolated_voltages(rMechanicsMesh.GetNumNodes());
    Vec voltage_coarse = NULL;

    for(unsigned time_step=0; time_step<num_timesteps; time_step++)
    {
        for (unsigned var_index = 0; var_index < rVariableNames.size(); var_index++)
        {
            std::string var_name = rVariableNames[var_index];
            // read
            reader.GetVariableOverNodes(voltage, var_name, time_step);
            ReplicatableVector voltage_repl(voltage);

            // interpolate
            for(unsigned i=0; i<mesh_pair.rGetElementsAndWeights().size(); i++)
            {
                double interpolated_voltage = 0;

                Element<DIM,DIM>& element = *(rElectricsMesh.GetElement(mesh_pair.rGetElementsAndWeights()[i].ElementNum));
                for(unsigned node_index = 0; node_index<element.GetNumNodes(); node_index++)
                {
                    unsigned global_node_index = element.GetNodeGlobalIndex(node_index);
                    interpolated_voltage += voltage_repl[global_node_index]*mesh_pair.rGetElementsAndWeights()[i].Weights(node_index);
                }

                interpolated_voltages[i] = interpolated_voltage;
            }

            if(voltage_coarse!=NULL)
            {
                PetscTools::Destroy(voltage_coarse);
            }
            voltage_coarse = PetscTools::CreateVec(interpolated_voltages);
            // write
            p_writer->PutVector(columns_id[var_index], voltage_coarse);
        }
        p_writer->PutUnlimitedVariable(time_step);
        p_writer->AdvanceAlongUnlimitedDimension();
    }

    if(voltage_coarse!=NULL)
    {
        PetscTools::Destroy(voltage);
        PetscTools::Destroy(voltage_coarse);
    }

    // delete to flush
    delete p_writer;

    // Convert the new data to CMGUI format.
    // alter the directory in HeartConfig as that is where Hdf5ToCmguiConverter decides
    // where to output
    std::string config_directory = HeartConfig::Instance()->GetOutputDirectory();
    HeartConfig::Instance()->SetOutputDirectory(directory);
    Hdf5ToCmguiConverter<DIM,DIM> converter(FileFinder(directory, RelativeTo::ChasteTestOutput),
                                            "voltage_mechanics_mesh",
                                            &rMechanicsMesh,
                                            false);
    HeartConfig::Instance()->SetOutputDirectory(config_directory);
}

    // 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
        }
    }