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

    void TestNodeExchange() throw(Exception)
    {
        HeartConfig::Instance()->Reset();

        HeartConfig::Instance()->Reset();
        DistributedTetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(0.1, 1.0); // [0,1] with h=0.1, ie 11 node mesh

        MyCardiacCellFactory cell_factory;
        cell_factory.SetMesh(&mesh);

        MonodomainTissue<1> monodomain_tissue( &cell_factory, true );


        if ( PetscTools::GetNumProcs() == 1 )
        {
            TS_ASSERT_EQUALS( mesh.GetNumHaloNodes(), 0u );
        }
        else
        {
            if ( PetscTools::AmMaster() || PetscTools::AmTopMost() )
            {
                TS_ASSERT_EQUALS( mesh.GetNumHaloNodes(), 1u );
            }
            else
            {
                TS_ASSERT_EQUALS( mesh.GetNumHaloNodes(), 2u );
            }
        }

        for (DistributedTetrahedralMesh<1,1>::HaloNodeIterator it=mesh.GetHaloNodeIteratorBegin();
                it != mesh.GetHaloNodeIteratorEnd();
                ++it)
        {
            AbstractCardiacCellInterface* cell = monodomain_tissue.GetCardiacCellOrHaloCell( (*it)->GetIndex() );
            TS_ASSERT_DELTA(cell->GetStimulus(0.001),0,1e-10);
        }

        if ( PetscTools::AmMaster() )
        {
            // Master owns node 0
            AbstractCardiacCellInterface* cell = monodomain_tissue.GetCardiacCellOrHaloCell(0);
            TS_ASSERT_DELTA(cell->GetStimulus(0.001), -80.0, 1e-10);
        }
        else
        {
            // Zero is not halo owned by any process (unless we have a lot of them).
            TS_ASSERT_THROWS_CONTAINS(monodomain_tissue.GetCardiacCellOrHaloCell(0),
                                      "Requested node/halo 0 does not belong to processor ");
        }
    }

    /*
     * This is the same as TestConductionVelocityConvergesFasterWithSvi1d with i=2, but solves in two parts.
     * If that test changes, check the hardcoded values here!
     */
    void TestArchiving()
    {
        std::string archive_dir = "monodomain_svi_archive";
        std::string archive_file = "monodomain_svi.arch";
        std::string output_dir = "monodomain_svi_output";

        { // Save
            HeartConfig::Instance()->SetSimulationDuration(0.1); //ms
            HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.01, 0.01);
            HeartConfig::Instance()->SetOutputDirectory(output_dir);
            HeartConfig::Instance()->SetOutputFilenamePrefix("results");
            HeartConfig::Instance()->SetUseStateVariableInterpolation();

            DistributedTetrahedralMesh<1,1> mesh;
            mesh.ConstructRegularSlabMesh(0.02, 1.0);

            BlockCellFactory<1> cell_factory;
            MonodomainProblem<1> monodomain_problem( &cell_factory );
            monodomain_problem.SetMesh(&mesh);
            monodomain_problem.Initialise();
            monodomain_problem.Solve();

            CardiacSimulationArchiver<MonodomainProblem<1> >::Save(monodomain_problem, archive_dir);
        }

        HeartConfig::Instance()->Reset();

        { // Load
            HeartConfig::Instance()->SetUseStateVariableInterpolation(false); // Just in case...
            MonodomainProblem<1> *p_monodomain_problem = CardiacSimulationArchiver<MonodomainProblem<1> >::Load(archive_dir);

            HeartConfig::Instance()->SetSimulationDuration(4.0); //ms
            p_monodomain_problem->Solve();

            ReplicatableVector final_voltage;
            final_voltage.ReplicatePetscVector(p_monodomain_problem->GetSolution());
            double probe_voltage = final_voltage[15];
            TS_ASSERT_DELTA(probe_voltage, 17.3131, 1e-3);

            delete p_monodomain_problem;
        }
    }

    void TestConductionVelocityConvergesFasterWithSvi1d()
    {
        double h[3] = {0.001,0.01,0.02};
        unsigned probe_node_index[3] = {300, 30, 15};
        unsigned number_of_nodes[3] = {1001, 101, 51};
        std::vector<double> conduction_vel_ici(3);
        std::vector<double> conduction_vel_svi(3);

        ReplicatableVector final_voltage_ici;
        ReplicatableVector final_voltage_svi;

        //HeartConfig::Instance()->SetUseRelativeTolerance(1e-8);
        HeartConfig::Instance()->SetSimulationDuration(4.0); //ms
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.01, 0.01);

        for (unsigned i=0; i<3; i++)
        {
            // ICI - ionic current interpolation - the default
            {
                DistributedTetrahedralMesh<1,1> mesh;
                mesh.ConstructRegularSlabMesh(h[i], 1.0);
                TS_ASSERT_EQUALS(mesh.GetNumNodes(), number_of_nodes[i]);

                //Double check (for later) that the indexing is as expected
                if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal( probe_node_index[i] ))
                {
                    TS_ASSERT_DELTA(mesh.GetNode( probe_node_index[i] )->rGetLocation()[0], 0.3, 1e-8);
                }
                std::stringstream output_dir;
                output_dir << "MonodomainIci_" << h[i];
                HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
                HeartConfig::Instance()->SetOutputFilenamePrefix("results");

                // need to have this for i=1,2 cases!!
                HeartConfig::Instance()->SetUseStateVariableInterpolation(false);

                BlockCellFactory<1> cell_factory;
                MonodomainProblem<1> monodomain_problem( &cell_factory );
                monodomain_problem.SetMesh(&mesh);
                monodomain_problem.Initialise();

                monodomain_problem.Solve();

                final_voltage_ici.ReplicatePetscVector(monodomain_problem.GetSolution());

//// see #1633
//// end time needs to be increased for these (say, to 7ms)
//                Hdf5DataReader simulation_data(OutputFileHandler::GetChasteTestOutputDirectory() + output_dir.str(),
//                                               "results", false);
//                PropagationPropertiesCalculator ppc(&simulation_data);
//                unsigned node_at_0_04 = (unsigned)round(0.04/h[i]);
//                unsigned node_at_0_40 = (unsigned)round(0.40/h[i]);
//                assert(fabs(mesh.GetNode(node_at_0_04)->rGetLocation()[0]-0.04)<1e-6);
//                assert(fabs(mesh.GetNode(node_at_0_40)->rGetLocation()[0]-0.40)<1e-6);
//                conduction_vel_ici[i] = ppc.CalculateConductionVelocity(node_at_0_04,node_at_0_40,0.36);
//                std::cout << "conduction_vel_ici = " << conduction_vel_ici[i] << "\n";
            }

            // SVI - state variable interpolation
            {
                DistributedTetrahedralMesh<1,1> mesh;
                mesh.ConstructRegularSlabMesh(h[i], 1.0);

                //Double check (for later) that the indexing is as expected
                if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal( probe_node_index[i] ))
                {
                    TS_ASSERT_DELTA(mesh.GetNode( probe_node_index[i] )->rGetLocation()[0], 0.3, 1e-8);
                }
                std::stringstream output_dir;
                output_dir << "MonodomainSvi_" << h[i];
                HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
                HeartConfig::Instance()->SetOutputFilenamePrefix("results");

                HeartConfig::Instance()->SetUseStateVariableInterpolation();

                BlockCellFactory<1> cell_factory;
                MonodomainProblem<1> monodomain_problem( &cell_factory );
                monodomain_problem.SetMesh(&mesh);
                monodomain_problem.Initialise();

                monodomain_problem.Solve();

                final_voltage_svi.ReplicatePetscVector(monodomain_problem.GetSolution());

//                Hdf5DataReader simulation_data(OutputFileHandler::GetChasteTestOutputDirectory() + output_dir.str(),
//                                               "results", false);
//                PropagationPropertiesCalculator ppc(&simulation_data);
//                unsigned node_at_0_04 = (unsigned)round(0.04/h[i]);
//                unsigned node_at_0_40 = (unsigned)round(0.40/h[i]);
//                assert(fabs(mesh.GetNode(node_at_0_04)->rGetLocation()[0]-0.04)<1e-6);
//                assert(fabs(mesh.GetNode(node_at_0_40)->rGetLocation()[0]-0.40)<1e-6);
//                conduction_vel_svi[i] = ppc.CalculateConductionVelocity(node_at_0_04,node_at_0_40,0.36);
//                std::cout << "conduction_vel_svi = " << conduction_vel_svi[i] << "\n";
            }

            if (i==0) // finest mesh
            {
                for (unsigned j=0; j<final_voltage_ici.GetSize(); j++)
                {
                    // visually checked they agree at this mesh resolution, and chosen tolerance from results
//.........这里部分代码省略.........

    void Test1dApd() throw(Exception)
    {
        HeartConfig::Instance()->SetPrintingTimeStep(1.0);
        HeartConfig::Instance()->SetSimulationDuration(400); //ms

        DistributedTetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(0.01, 1.0); // h=0.01cm, width=1cm

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory(-600.0*1000);

        //////////////////////////////////////////////////////////////////////////
        // run original simulation - no adaptivity, dt=0.01 all the way through
        //////////////////////////////////////////////////////////////////////////
        HeartConfig::Instance()->SetOutputDirectory("MonoWithTimeAdaptivity1dLong/OrigNoAdapt");
        MonodomainProblem<1> problem(&cell_factory);
        problem.SetMesh(&mesh);

        problem.Initialise();
        problem.Solve();

        HeartEventHandler::Headings();
        HeartEventHandler::Report();

        //////////////////////////////////////////////////////////////////////////
        // run adaptive simulation - dt=0.01 for first 2ms, then dt=1
        //////////////////////////////////////////////////////////////////////////
        HeartConfig::Instance()->SetOutputDirectory("MonoWithTimeAdaptivity1dLong/SimpleAdapt");
        MonodomainProblem<1> adaptive_problem(&cell_factory);
        adaptive_problem.SetMesh(&mesh);

        FixedTimeAdaptivityController controller(25);
        adaptive_problem.SetUseTimeAdaptivityController(true, &controller);
        adaptive_problem.Initialise();
        adaptive_problem.Solve();

        HeartEventHandler::Headings();
        HeartEventHandler::Report();

        Hdf5DataReader reader_no_adapt("MonoWithTimeAdaptivity1dLong/OrigNoAdapt","SimulationResults");
        Hdf5DataReader reader_adapt("MonoWithTimeAdaptivity1dLong/SimpleAdapt","SimulationResults");

        unsigned num_timesteps = reader_no_adapt.GetUnlimitedDimensionValues().size();
        assert(num_timesteps == reader_adapt.GetUnlimitedDimensionValues().size());

        DistributedVectorFactory factory(mesh.GetNumNodes());
        Vec voltage_no_adapt = factory.CreateVec();
        Vec voltage_adapt = factory.CreateVec();

        Vec difference;
        VecDuplicate(voltage_adapt, &difference);

        for (unsigned timestep=0; timestep<num_timesteps; timestep++)
        {
            reader_no_adapt.GetVariableOverNodes(voltage_no_adapt, "V", timestep);
            reader_adapt.GetVariableOverNodes(voltage_adapt, "V", timestep);

            PetscVecTools::WAXPY(difference, -1.0, voltage_adapt, voltage_no_adapt);
            double l_inf_norm;
            VecNorm(difference, NORM_INFINITY, &l_inf_norm);

            //std::cout << l_inf_norm << "\n";
            if (timestep < 25)
            {
                TS_ASSERT_DELTA(l_inf_norm, 0.0, 1e-10); // first 25 ms, there should be no difference
            }
            else
            {
                TS_ASSERT_DELTA(l_inf_norm, 0.0, 2.25); // the difference is at most ~2mv, which occurs during the downstroke
            }
        }

        PetscTools::Destroy(voltage_no_adapt);
        PetscTools::Destroy(voltage_adapt);
    }

    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 TestSimpleSimulation() throw(Exception)
    {
        /*Simulation parameters*/
        HeartConfig::Instance()->SetSimulationDuration(0.7); //ms (falls over after this)
        HeartConfig::Instance()->SetUseAbsoluteTolerance(1e-6);
        //HeartConfig::Instance()->SetOdeTimeStep(0.01);

        const double width = 0.1;
        const double height = 0.1;
        const double depth = 0.1;

        const unsigned num_elem_x = 8;
        const double space_step = width/num_elem_x;

        /* Make the mesh*/
        DistributedTetrahedralMesh<3,3> mesh;
        mesh.ConstructRegularSlabMesh(space_step, width, height, depth);

        /*Create a cell factory of the type we defined above. */
        GeneralPlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 3> cell_factory(num_elem_x, width);

        /* monodomain problem class using (a pointer to) the cell factory */
        BidomainProblem<3> problem( &cell_factory );
        problem.SetMesh(&mesh);

        /*
        * HOW_TO_TAG Cardiac/Problem definition
        * Set discrete '''cuboid''' areas to have heterogeneous (intra- and/or extra-cellular) conductivity tensors.
        */
        std::vector<ChasteCuboid<3> > input_areas;
        std::vector< c_vector<double,3> > intra_conductivities;
        std::vector< c_vector<double,3> > extra_conductivities;
        ChastePoint<3> corner_a(width/2, 0, 0);
        ChastePoint<3> corner_b(width, height, depth);

        input_areas.push_back(ChasteCuboid<3> (corner_a, corner_b));
        //within the cuboid
        intra_conductivities.push_back( Create_c_vector(0.1, 0.1, 0.1) );
        extra_conductivities.push_back( Create_c_vector(0.0, 0.0, 0.0) );
        //This test should *fail* if you comment out the following line
        //(which blocks conductivity on the RHS of the slab).
        HeartConfig::Instance()->SetConductivityHeterogeneities(input_areas, intra_conductivities, extra_conductivities);

        //elsewhere
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(1.2, 1.2, 1.2));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(1.2, 1.2, 1.2));

        /* set  parameters*/
        // HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1.0);
        // HeartConfig::Instance()->SetCapacitance(1.0);

         /* Output Directory and prefix (for the hdf5 file), relative to CHASTE_TEST_OUTPUT*/
        HeartConfig::Instance()->SetOutputDirectory("slab_results_het_halfcond");
        HeartConfig::Instance()->SetOutputFilenamePrefix("Slab_small");

        /* Initialise the problem*/
        problem.Initialise();

        /* Solve the PDE monodomain equaion*/
        problem.Solve();

        ReplicatableVector voltage_replicated(problem.GetSolution());
        TS_ASSERT_EQUALS(mesh.GetNumNodes() * 2, voltage_replicated.GetSize());
        unsigned lo, hi;
        lo = mesh.GetDistributedVectorFactory()->GetLow();
        hi = mesh.GetDistributedVectorFactory()->GetHigh();

        for (unsigned i=lo; i<hi; i++)
        {
            double x = mesh.GetNode(i)->rGetLocation()[0];
            if (x<width/2)
            {
                 //Left side is stimulated
                 TS_ASSERT_LESS_THAN(-71.0,voltage_replicated[2 * i]);
            }
            else if (x>width/2)
            {
                //Right side is blocked
                TS_ASSERT_LESS_THAN(voltage_replicated[2 * i],-82.0);
            }
        }
    }

    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 TestSolveCellSystemsInclUpdateVoltageWithNodeExchange() throw(Exception)
    {
        HeartConfig::Instance()->Reset();
        DistributedTetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(1.0, 1.0); // [0,1] with h=0.1, ie 11 node mesh

        MyCardiacCellFactory cell_factory;
        cell_factory.SetMesh(&mesh);

        MonodomainTissue<1> monodomain_tissue( &cell_factory, true );

        Vec voltage = PetscTools::CreateAndSetVec(2, -81.4354); // something that isn't resting potential
        monodomain_tissue.SolveCellSystems(voltage, 0, 1, false); // solve for 1ms without updating the voltage

        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(0)) // Is node 0 locally owned?
        {
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCell(0)->GetVoltage(), -81.4354, 1e-3);
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCellOrHaloCell(1)->GetVoltage(), -81.4354, 1e-3);
        }

        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(1)) // Is node 1 locally owned?
        {
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCellOrHaloCell(0)->GetVoltage(), -81.4354, 1e-3);
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCell(1)->GetVoltage(), -81.4354, 1e-3);
        }

        Vec voltage2 = PetscTools::CreateAndSetVec(2, -75);
        monodomain_tissue.SolveCellSystems(voltage2, 1, 2, true); // solve another ms, using this new voltage, but now updating the voltage too
        ReplicatableVector voltage2_repl(voltage2); // should have changed following solve

        // check the new voltage in the cell is NEAR -75 (otherwise the passed in voltage wasn't used, but
        // NOT EXACTLY -75, ie that the voltage was solved for.
        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(0)) // Is node 0 locally owned?
        {
            // check has been updated
            TS_ASSERT_DIFFERS(monodomain_tissue.GetCardiacCell(0)->GetVoltage(), -75);
            // check near -75
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCell(0)->GetVoltage(), -75, 2.0); // within 2mV
            // check the passed in voltage was updated
            TS_ASSERT_DELTA(voltage2_repl[0], monodomain_tissue.GetCardiacCell(0)->GetVoltage(), 1e-10);
        }
        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(1)) // Is node 1 locally owned?
        {
            TS_ASSERT_DIFFERS(monodomain_tissue.GetCardiacCell(1)->GetVoltage(), -75);
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCell(1)->GetVoltage(), -75, 2.0); // within 2mV
            TS_ASSERT_DELTA(voltage2_repl[1], monodomain_tissue.GetCardiacCell(1)->GetVoltage(), 1e-10);
        }

        // now check the new voltages have been communicated
        // check the new voltage in the cell is NEAR -75 (otherwise the passed in voltage wasn't used, but
        // NOT EXACTLY -75, ie that the voltage was solved for.
        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(0)) // Is node 0 locally owned?
        {
            TS_ASSERT_DIFFERS(monodomain_tissue.GetCardiacCellOrHaloCell(1)->GetVoltage(), -75);
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCellOrHaloCell(1)->GetVoltage(), -75, 2.0); // within 2mV
            TS_ASSERT_DELTA(voltage2_repl[1], monodomain_tissue.GetCardiacCellOrHaloCell(1)->GetVoltage(), 1e-10);
        }
        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(1)) // Is node 1 locally owned?
        {
            TS_ASSERT_DIFFERS(monodomain_tissue.GetCardiacCellOrHaloCell(0)->GetVoltage(), -75);
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCellOrHaloCell(0)->GetVoltage(), -75, 2.0); // within 2mV
            TS_ASSERT_DELTA(voltage2_repl[0], monodomain_tissue.GetCardiacCellOrHaloCell(0)->GetVoltage(), 1e-10);
        }

        PetscTools::Destroy(voltage);
        PetscTools::Destroy(voltage2);
    }

    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

    }

    void TestConductionVelocityConvergesFasterWithSvi1d() throw(Exception)
    {
        double h[3] = {0.001,0.01,0.02};
        std::vector<double> conduction_vel_nci(3);
        std::vector<double> conduction_vel_svi(3);

        ReplicatableVector final_solution_ici;
        ReplicatableVector final_solution_svi;
        ReplicatableVector final_solution_svit;

        HeartConfig::Instance()->SetSimulationDuration(4.0); //ms
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.01, 0.01, 0.01);

        for(unsigned i=0; i<3; i++)
        {
            // ICI - ionic current interpolation - the default
            {
                TetrahedralMesh<1,1> mesh;
                mesh.ConstructRegularSlabMesh(h[i], 1.0);

                std::stringstream output_dir;
                output_dir << "BidomainIci_" << h[i];
                HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
                HeartConfig::Instance()->SetOutputFilenamePrefix("results");

                // need to have this for i=1,2 cases!!
                HeartConfig::Instance()->SetUseStateVariableInterpolation(false);

                BlockCellFactory<1> cell_factory;
                BidomainProblem<1> 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
            {
                DistributedTetrahedralMesh<1,1> mesh;
                mesh.ConstructRegularSlabMesh(h[i], 1.0);

                std::stringstream output_dir;
                output_dir << "BidomainSvi_" << h[i];
                HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
                HeartConfig::Instance()->SetOutputFilenamePrefix("results");

                HeartConfig::Instance()->SetUseStateVariableInterpolation();

                BlockCellFactory<1> cell_factory;
                BidomainProblem<1> bidomain_problem( &cell_factory );
                bidomain_problem.SetMesh(&mesh);
                bidomain_problem.Initialise();

                bidomain_problem.Solve();

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

            // SVIT - state variable interpolation on straight (not distributed) tetrahedral mesh
            {
                TetrahedralMesh<1,1> mesh;
                mesh.ConstructRegularSlabMesh(h[i], 1.0);

                std::stringstream output_dir;
                output_dir << "BidomainSviTet_" << h[i];
                HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
                HeartConfig::Instance()->SetOutputFilenamePrefix("results");

                HeartConfig::Instance()->SetUseStateVariableInterpolation();

                BlockCellFactory<1> cell_factory;
                BidomainProblem<1> bidomain_problem( &cell_factory );
                bidomain_problem.SetMesh(&mesh);
                bidomain_problem.Initialise();

                bidomain_problem.Solve();

                final_solution_svit.ReplicatePetscVector(bidomain_problem.GetSolution());
            }
            double voltage_at_0_03_finest_mesh;
            if(i==0) // finest mesh
            {
                for(unsigned j=0; j<final_solution_ici.GetSize(); j++)
                {
                    // visually checked they agree at this mesh resolution, and chosen tolerance from results
                    TS_ASSERT_DELTA(final_solution_ici[j], final_solution_svi[j], 0.35);
                    TS_ASSERT_DELTA(final_solution_svit[j], final_solution_svi[j], 1e-8);

                    if(j%2==0 /* just look at voltage */ && final_solution_ici[j]>-80)
                    {
                        // shouldn't be exactly equal, as long as away from resting potential
                        TS_ASSERT_DIFFERS(final_solution_ici[j], final_solution_svi[j]);
                    }
                }

                voltage_at_0_03_finest_mesh = final_solution_ici[600];
                TS_ASSERT_DELTA(voltage_at_0_03_finest_mesh, -65.2218, 1e-2); //hardcoded value
            }
//.........这里部分代码省略.........