C++ DistributedTetrahedralMesh::GetElementIteratorEnd方法代码示例

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

        c_matrix<double, 3, 3> dummy;
        double jacobian_det = 0.0;

        double scaled_volume = 0.0;
        for (AbstractTetrahedralMesh<3, 3>::ElementIterator iter = mesh.GetElementIteratorBegin();
                iter != mesh.GetElementIteratorEnd();
                ++iter)
        {
            iter->CalculateJacobian(dummy, jacobian_det);
            scaled_volume += jacobian_det;
        }
        TS_ASSERT_DELTA(scaled_volume, 1.0*6, 1e-6);

        mesh.RotateX(M_PI);
        //mesh.Translate(100.0, 0.0, 0.0);

        double scaled_volume_after = 0.0;
        for (AbstractTetrahedralMesh<3, 3>::ElementIterator iter = mesh.GetElementIteratorBegin();
                iter != mesh.GetElementIteratorEnd();
                ++iter)
        {
            iter->CalculateJacobian(dummy, jacobian_det);
            scaled_volume_after += jacobian_det;
        }
        TS_ASSERT_DELTA(scaled_volume_after, 1.0*6, 1e-6);

    }

    void TestConductivityModifier() throw(Exception)
    {
        /*
         * Generate a mesh.
         */
        DistributedTetrahedralMesh<2,2> mesh;
        mesh.ConstructRegularSlabMesh(0.5, 1.0, 0.5); // Mesh has 4 elements

        /*
         * Here we're using a trivial cell factory for simplicity, but usually you'll provide your own one.
         * Set up the problem with the factory as usual.
         */
        ZeroStimulusCellFactory<CellLuoRudy1991FromCellML,2> cell_factory;
        BidomainProblem<2> bidomain_problem( &cell_factory );
        bidomain_problem.SetMesh( &mesh );

        /*
         * We need to apply the modifier directly to the tissue, which comes from the problem, but is only
         * accessible after `Initialise()`, so let's do that now.
         */
        bidomain_problem.Initialise();
        BidomainTissue<2>* p_bidomain_tissue = bidomain_problem.GetBidomainTissue();

        /*
         * Get the original conductivity tensor values. We haven't set them using
         * `HeartConfig->SetIntra/ExtracellularConductivities` so they'll just be the defaults.
         *
         * The first argument below is the element ID (we just check the first element we own here). The second accesses
         * the diagonal elements. We just do (0,0), as (1,1) should be the same (no fibre orientation).
         * Off-diagonal elements will be 0.
         *
         * As we don't have many elements, when we run on more than two processors some processors
         * will not own any elements. We only try to access the conductivity tensors if the process
         * owns at least one element.
         *
         * We then check that we have the correct (default) conductivity values.
         */
        double orig_intra_conductivity_0 = 0.0;
        double orig_extra_conductivity_0 = 0.0;
        if (mesh.GetElementIteratorBegin() != mesh.GetElementIteratorEnd())
        {
            unsigned first_element = mesh.GetElementIteratorBegin()->GetIndex();
            orig_intra_conductivity_0 = p_bidomain_tissue->rGetIntracellularConductivityTensor(first_element)(0,0);
            orig_extra_conductivity_0 = p_bidomain_tissue->rGetExtracellularConductivityTensor(first_element)(0,0);

            TS_ASSERT_DELTA(orig_intra_conductivity_0, 1.75, 1e-9); // hard-coded using default
            TS_ASSERT_DELTA(orig_extra_conductivity_0, 7.0, 1e-9); // hard-coded using default
        }

        /*
         * Now we can make the modifier and apply it to the tissue using `SetConductivityModifier`.
         */
        SimpleConductivityModifier modifier;
        p_bidomain_tissue->SetConductivityModifier( &modifier );

        /*
         * To confirm that the conductivities have changed, let's iterate over all elements owned by this process
         * and check their conductivity against what we expect.
         */
        for (AbstractTetrahedralMesh<2,2>::ElementIterator elt_iter=mesh.GetElementIteratorBegin();
             elt_iter!=mesh.GetElementIteratorEnd();
             ++elt_iter)
        {
            unsigned index = elt_iter->GetIndex();
            if (index == 0u)
            {
                TS_ASSERT_DELTA(p_bidomain_tissue->rGetIntracellularConductivityTensor(0)(0,0), 3.14, 1e-9);
                TS_ASSERT_DELTA(p_bidomain_tissue->rGetExtracellularConductivityTensor(0)(0,0), 3.14, 1e-9);
                TS_ASSERT_DELTA(p_bidomain_tissue->rGetIntracellularConductivityTensor(0)(1,1), 0.707, 1e-9);
                TS_ASSERT_DELTA(p_bidomain_tissue->rGetExtracellularConductivityTensor(0)(1,1), 0.707, 1e-9);
            }
            else
            {
                TS_ASSERT_DELTA(p_bidomain_tissue->rGetIntracellularConductivityTensor(index)(0,0), 1.0*index*orig_intra_conductivity_0, 1e-9);
                TS_ASSERT_DELTA(p_bidomain_tissue->rGetExtracellularConductivityTensor(index)(0,0), 1.5*index*orig_extra_conductivity_0, 1e-9);
            }
        }
    }

    void TestWithBathAndElectrodes()
    {
        /* First, set the end time and output info. In this simulation
         * we'll explicitly read the mesh, alter it, then pass it
         * to the problem class, so we don't set the mesh file name.
         */
        HeartConfig::Instance()->SetSimulationDuration(3.0);  //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainTutorialWithBath");
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");

        /* Bath problems seem to require decreased ODE timesteps.
         */
        HeartConfig::Instance()->SetOdeTimeStep(0.001);  //ms

        /* Use the {{{PlaneStimulusCellFactory}}} to define a set
         * of Luo-Rudy cells. We pass the stimulus magnitude as 0.0
         * as we don't want any stimulated cells.
         */
        PlaneStimulusCellFactory<CellLuoRudy1991FromCellMLBackwardEuler,2> cell_factory(0.0);

        /*
         * Now, we load up a rectangular mesh (in triangle/tetgen format), done as follows,
         * using {{{TrianglesMeshReader}}}.  Note that we use a distributed mesh, so the data
         * is shared among processes if run in parallel.
         */
        TrianglesMeshReader<2,2> reader("mesh/test/data/2D_0_to_1mm_400_elements");
        DistributedTetrahedralMesh<2,2> mesh;
        mesh.ConstructFromMeshReader(reader);

        /*
         * In most simulations there is one valid tissue identifier and one valid bath identifier
         * (for elements).
         * One of these can be assigned to an element with
         *  * {{{mesh.GetElement(i)->SetAttribute(HeartRegionCode::GetValidTissueId());}}}
         *  * {{{mesh.GetElement(i)->SetAttribute(HeartRegionCode::GetValidBathId());}}}
         *
         * If we want heterogeneous conductivities outside the heart (for example for torso and blood)
         * then we will need different identifiers:
         */
        std::set<unsigned> tissue_ids;
        static unsigned tissue_id=0;
        tissue_ids.insert(tissue_id);

        std::set<unsigned> bath_ids;
        static unsigned bath_id1=1;
        bath_ids.insert(bath_id1);
        static unsigned bath_id2=2;
        bath_ids.insert(bath_id2);

        HeartConfig::Instance()->SetTissueAndBathIdentifiers(tissue_ids, bath_ids);

        /* In bath problems, each element has an attribute which must be set
         * to 0 (cardiac tissue) or 1 (bath). This can be done by having an
         * extra column in the element file (see the file formats documentation,
         * or for example
         * mesh/test/data/1D_0_to_1_10_elements_with_two_attributes.ele,
         * and note that the header in this file has 1 at the end to indicate that
         * the file defines an attribute for each element). We have read in a mesh
         * without this type of information set up, so we set it up manually,
         * by looping over elements and setting those more than 2mm from the centre
         * as bath elements (by default, the others are cardiac elements).
         */
        for (AbstractTetrahedralMesh<2,2>::ElementIterator iter = mesh.GetElementIteratorBegin();
             iter != mesh.GetElementIteratorEnd();
             ++iter)
        {
            double x = iter->CalculateCentroid()[0];
            double y = iter->CalculateCentroid()[1];
            if (sqrt((x-0.05)*(x-0.05) + (y-0.05)*(y-0.05)) > 0.02)
            {
                if (y<0.05)
                {
                    //Outside circle on the bottom
                    iter->SetAttribute(bath_id1);
                }
                else
                {
                    //Outside circle on the top
                    iter->SetAttribute(bath_id2);
                }
            }
            else
            {
                //IDs default to 0, but we want to be safe
                iter->SetAttribute(tissue_id);
            }
        }

        /* HOW_TO_TAG Cardiac/Problem definition
         * Tell Chaste that a mesh has been modified
         *
         * Since we have modified the mesh by setting element attributes, we need to inform Chaste of this fact.
         * If we do not, problems will arise when [wiki:UserTutorials/CardiacCheckpointingAndRestarting checkpointing],
         * since the code that saves the simulation state will assume that it can just reuse the original mesh files,
         * and thus won't save the new element attributes.
         *
         * (Some mesh modifications, that use methods on the mesh class directly, will automatically record that
         * the mesh has been modified.  Since we're just modifying elements, this information isn't propagated at
         * present.)
         */
//.........这里部分代码省略.........

    void Test2dBathMultipleBathConductivities() throw (Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(2.0);  //ms
        HeartConfig::Instance()->SetOutputDirectory("BidomainBath2dMultipleBathConductivities");
        HeartConfig::Instance()->SetOutputFilenamePrefix("bidomain_bath_2d");

        HeartConfig::Instance()->SetOdeTimeStep(0.001);  //ms ???

        std::set<unsigned> tissue_ids;
        tissue_ids.insert(0); // Same as default value defined in HeartConfig

        std::set<unsigned> bath_ids;
        bath_ids.insert(2);
        bath_ids.insert(3);
        bath_ids.insert(4);
        HeartConfig::Instance()->SetTissueAndBathIdentifiers(tissue_ids, bath_ids);

        // need to create a cell factory but don't want any intra stim, so magnitude
        // of stim is zero.
        c_vector<double,2> centre;
        centre(0) = 0.05; // cm
        centre(1) = 0.05; // cm
        BathCellFactory<2> cell_factory( 0.0, centre);

        BidomainWithBathProblem<2> bidomain_problem( &cell_factory );

        DistributedTetrahedralMesh<2,2> mesh;

        mesh.ConstructRegularSlabMesh(0.05, 0.9, 0.9);

        // set the x<0.25 and x>0.75 regions as the bath region
        for (AbstractTetrahedralMesh<2,2>::ElementIterator iter = mesh.GetElementIteratorBegin();
             iter != mesh.GetElementIteratorEnd();
             ++iter)
        {
            double x = iter->CalculateCentroid()[0];
            double y = iter->CalculateCentroid()[1];
            if( (x>0.3) && (x<0.6) && (y>0.3) && (y<0.6) )
            {
                iter->SetAttribute(0);
            }
            else
            {
                if (y<0.2)
                {
                    iter->SetAttribute(2);
                }
                else if (y<0.7)
                {
                    iter->SetAttribute(3);
                }
                else if (y<0.9)
                {
                    iter->SetAttribute(4);
                }
            }
        }

        std::map<unsigned, double> multiple_bath_conductivities;
        multiple_bath_conductivities[2] = 7.0;
        multiple_bath_conductivities[3] = 1.0;
        multiple_bath_conductivities[4] = 0.001;

        HeartConfig::Instance()->SetBathMultipleConductivities(multiple_bath_conductivities);

        double boundary_flux = -3.0e3;
        double start_time = 0.0;
        double duration = 1.0; // of the stimulus, in ms

        HeartConfig::Instance()->SetElectrodeParameters(false, 0, boundary_flux, start_time, duration);

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

        bidomain_problem.Solve();

        DistributedVector distributed_solution = bidomain_problem.GetSolutionDistributedVector();
        DistributedVector::Stripe voltage(distributed_solution, 0);

        /*
         * We are checking the last time step. This test will only make sure that an AP is triggered.
         */
        bool ap_triggered = false;

        for (DistributedVector::Iterator index = distributed_solution.Begin();
             index!= distributed_solution.End();
             ++index)
        {
            // test V = 0 for all bath nodes and that an AP is triggered in the tissue
            if (HeartRegionCode::IsRegionBath( mesh.GetNode(index.Global)->GetRegion() )) // bath
            {
                TS_ASSERT_DELTA(voltage[index], 0.0, 1e-12);
            }
            else if (voltage[index] > 0.0)//at the last time step
            {
                ap_triggered = true;
            }
        }

        TS_ASSERT(PetscTools::ReplicateBool(ap_triggered));
//.........这里部分代码省略.........

    void TestBidomainWithBathWithSvi() throw(Exception)
    {
        /* Make a 4x4 node mesh and set two interior elements to be bath elements */
        DistributedTetrahedralMesh<2,2> mesh;
        mesh.ConstructRegularSlabMesh(0.04, 0.12, 0.12);

        for (AbstractTetrahedralMesh<2,2>::ElementIterator iter=mesh.GetElementIteratorBegin();
             iter != mesh.GetElementIteratorEnd();
             ++iter)
        {
            unsigned element_index = iter->GetIndex();

            if ( element_index==10 || element_index==17 )
            {
                iter->SetAttribute(HeartRegionCode::GetValidBathId());
            }
            else
            {
                iter->SetAttribute(HeartRegionCode::GetValidTissueId());
            }
        }

        HeartConfig::Instance()->SetSimulationDuration(10.0); // ms
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.001, 0.025, 0.25);
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(1.75, 0.17));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(7.0, 0.7));

        ReplicatableVector final_solution_ici;
        ReplicatableVector final_solution_svi;

        // ICI - ionic current interpolation (the default)
        {
            HeartConfig::Instance()->SetOutputDirectory("BidomainWithBathIci2d");
            HeartConfig::Instance()->SetOutputFilenamePrefix("results");
            HeartConfig::Instance()->SetUseStateVariableInterpolation(false);

            BathCellFactory cell_factory;
            BidomainWithBathProblem<2> bidomain_problem( &cell_factory );
            bidomain_problem.SetMesh(&mesh);
            bidomain_problem.Initialise();
            bidomain_problem.Solve();

            final_solution_ici.ReplicatePetscVector(bidomain_problem.GetSolution());
        }

        // SVI - state variable interpolation
        {
            HeartConfig::Instance()->SetOutputDirectory("BidomainWithBathSvi2d");
            HeartConfig::Instance()->SetOutputFilenamePrefix("results");
            HeartConfig::Instance()->SetUseStateVariableInterpolation(true);

            BathCellFactory cell_factory;
            BidomainWithBathProblem<2> bidomain_problem( &cell_factory );
            bidomain_problem.SetMesh(&mesh);
            bidomain_problem.Initialise();
            bidomain_problem.Solve();

            final_solution_svi.ReplicatePetscVector(bidomain_problem.GetSolution());
        }

        // ICI
        TS_ASSERT_DELTA(final_solution_ici[15*2], 7.0918, 2e-3); // Node 15 phi_i
        TS_ASSERT_DELTA(final_solution_ici[15*2+1], 0.0401, 1e-3); // Node 15 phi_e
        // SVI
        TS_ASSERT_DELTA(final_solution_svi[15*2], 10.6217, 2e-3); // Node 15 phi_i
        TS_ASSERT_DELTA(final_solution_svi[15*2+1], -0.0180, 1e-3); // Node 15 phi_e

    }