C++ DistributedTetrahedralMesh类代码示例

    void TestConstructStreeterOnRightWedge() throw(Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/human_wedge_mesh/HumanWedgeMesh");
        std::string epi_face_file = "heart/test/data/human_wedge_mesh/epi.tri";
        std::string endo_face_file = "heart/test/data/human_wedge_mesh/endo.tri";

        DistributedTetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        StreeterFibreGenerator<3> fibre_generator(mesh);

        //Assume we are in the left ventricle
        fibre_generator.SetSurfaceFiles(epi_face_file, endo_face_file, "", true);

        fibre_generator.SetApexToBase(0);
        fibre_generator.SetWriteFileAsBinary();

        OutputFileHandler handler("human_wedge_mesh", false);

        fibre_generator.WriteData(handler, "HumanWedgeMeshRight.ortho");

        FileFinder fibre_file1 = handler.FindFile("HumanWedgeMeshRight.ortho");
        FileFinder fibre_file2("heart/test/data/human_wedge_mesh/HumanWedgeMeshRight.ortho", RelativeTo::ChasteSourceRoot);

        CompareGeneratedWithReferenceFile(fibre_file1, ORTHO, fibre_file2, ORTHO);
    }

    void TestSetLogInfo() throw (Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/box_shaped_heart/box_heart");
        std::string epi_face_file = "heart/test/data/box_shaped_heart/epi.tri";
        std::string rv_face_file = "heart/test/data/box_shaped_heart/rv.tri";
        std::string lv_face_file = "heart/test/data/box_shaped_heart/lv.tri";

        DistributedTetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        StreeterFibreGenerator<3> fibre_generator(mesh);
        fibre_generator.SetSurfaceFiles(epi_face_file, rv_face_file, lv_face_file, false);
        fibre_generator.SetApexToBase(0);

        OutputFileHandler handler("shorter_streeter_loginfo");

        fibre_generator.WriteData(handler, "box_heart.ortho");

        FileFinder node_regions_file = handler.FindFile("node_regions.data");
        FileFinder wall_thickness_file = handler.FindFile("wall_thickness.data");
        FileFinder averaged_thickness_file = handler.FindFile("averaged_thickness.data");

        TS_ASSERT_EQUALS(node_regions_file.IsFile(), false);
        TS_ASSERT_EQUALS(wall_thickness_file.IsFile(), false);
        TS_ASSERT_EQUALS(averaged_thickness_file.IsFile(), false);

        fibre_generator.SetLogInfo(true);

        fibre_generator.WriteData(handler, "box_heart.ortho");

        TS_ASSERT_EQUALS(node_regions_file.IsFile(), true);
        TS_ASSERT_EQUALS(wall_thickness_file.IsFile(), true);
        TS_ASSERT_EQUALS(averaged_thickness_file.IsFile(), true);
    }

    void TestSimpleOrthotropic() throw (Exception)
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/box_shaped_heart/box_heart");
        std::string epi_face_file = "heart/test/data/box_shaped_heart/epi.tri";
        std::string rv_face_file = "heart/test/data/box_shaped_heart/rv.tri";
        std::string lv_face_file = "heart/test/data/box_shaped_heart/lv.tri";

        DistributedTetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        StreeterFibreGenerator<3> fibre_generator(mesh);
        fibre_generator.SetSurfaceFiles(epi_face_file, rv_face_file, lv_face_file, false);
        fibre_generator.SetApexToBase(0);

        OutputFileHandler handler("shorter_streeter", false);
        fibre_generator.WriteData(handler, "box_heart.ortho");

        FileFinder fibre_file_ascii = handler.FindFile("box_heart.ortho");
        FileFinder fibre_file_reference("heart/test/data/box_shaped_heart/box_heart.ortho", RelativeTo::ChasteSourceRoot);

        CompareGeneratedWithReferenceFile(fibre_file_ascii, ORTHO, fibre_file_reference, ORTHO);

        fibre_generator.SetWriteFileAsBinary();
        fibre_generator.WriteData(handler, "box_heart_binary.ortho");

        FileFinder fibre_file_binary = handler.FindFile("box_heart_binary.ortho");

        CompareGeneratedWithReferenceFile(fibre_file_binary, ORTHO, fibre_file_reference, ORTHO);
    }

    void TestCheckForBathElementsNoDeadlock() throw (Exception)
    {
        HeartConfig::Instance()->SetSimulationDuration(1.0);  //ms
        HeartConfig::Instance()->SetOutputDirectory("bidomain_bath");
        HeartConfig::Instance()->SetOutputFilenamePrefix("BidomainLR91_1d");

        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> bidomain_cell_factory;
        BidomainWithBathProblem<1> bidomain_problem( &bidomain_cell_factory );

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

        try
        {
            mesh.GetElement(0)->SetAttribute(HeartRegionCode::GetValidBathId());
        }
        catch(Exception&)
        {
            // I don't own element 0
        }

        bidomain_problem.SetMesh(&mesh);

        // Fails because no bath
        TS_ASSERT_THROWS_NOTHING(bidomain_problem.Initialise());

        // Prevent an EventHandling exception in later tests
        HeartEventHandler::EndEvent(HeartEventHandler::EVERYTHING);
    }

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

    }

    /*
     * 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 TestExceptions()
    {
        TrianglesMeshReader<3,3> mesh_reader("heart/test/data/box_shaped_heart/box_heart");

        DistributedTetrahedralMesh<3,3> mesh;
        mesh.ConstructFromMeshReader(mesh_reader);

        StreeterFibreGenerator<3> fibre_generator(mesh);

        // No surfaces defined
        OutputFileHandler handler("streeter", false);
        TS_ASSERT_THROWS_THIS(fibre_generator.WriteData(handler, "file.fibres"),
                "Files defining the heart surfaces not set");

        // Wrong surface filename
        TS_ASSERT_THROWS_THIS(fibre_generator.SetSurfaceFiles("wrong_name", "wrong_name", "wrong_name", false),
                "Wrong surface definition file name wrong_name");

        std::string epi_face_file = "heart/test/data/box_shaped_heart/epi.tri";
        std::string rv_face_file = "heart/test/data/box_shaped_heart/rv.tri";
        std::string lv_face_file = "heart/test/data/box_shaped_heart/lv.tri";


        fibre_generator.SetSurfaceFiles(epi_face_file, rv_face_file, lv_face_file, false);
        OutputFileHandler shorter_handler("shorter_streeter", false);
        TS_ASSERT_THROWS_THIS(fibre_generator.WriteData(shorter_handler, "vector_not_set.ortho"),
            "Apex to base vector has not been set");
        TS_ASSERT_THROWS_THIS(fibre_generator.SetApexToBase(999),
            "Apex to base coordinate axis was out of range");

        c_vector<double, 3> axis;
        axis[0] = 0.0;
        axis[1] = 0.0;
        axis[2] = 0.0;
        TS_ASSERT_THROWS_THIS(fibre_generator.SetApexToBase(axis),
            "Apex to base vector should be non-zero");
        axis[1] = 42.0; //Will be normalised
        fibre_generator.SetApexToBase(axis);
    }

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

    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
            }
//.........这里部分代码省略.........

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

    /**
     * Tests archiving of the tissue object.
     * It creates one, changes the default values of some member variables and saves.
     * Then it tries to load from the archive and checks that the member variables are with the right values.
     */
    void TestSaveAndLoadExtendedBidomainTissue() throw (Exception)
    {
        HeartConfig::Instance()->Reset();
        // Archive settings
        FileFinder archive_dir("extended_tissue_archive", RelativeTo::ChasteTestOutput);
        std::string archive_file = "extended_bidomain_tissue.arch";

        bool cache_replication_saved = false;
        double saved_printing_timestep = 2.0;
        double default_printing_timestep = HeartConfig::Instance()->GetPrintingTimeStep();

        c_matrix<double, 3, 3> intra_tensor_before_archiving;
        c_matrix<double, 3, 3> intra_tensor_second_cell_before_archiving;
        c_matrix<double, 3, 3> extra_tensor_before_archiving;

        //creation and save
        {
            // This call is required to set the appropriate conductivity media and to make sure that HeartConfig
            // knows the mesh filename despite we use our own mesh reader.
            HeartConfig::Instance()->SetMeshFileName("mesh/test/data/cube_136_elements");

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

            UnStimulatedCellFactory first_cell;
            StimulatedCellFactory second_cell;
            ExtracellularStimulusFactory extra_factory;
            first_cell.SetMesh(&mesh);
            second_cell.SetMesh(&mesh);
            extra_factory.SetMesh(&mesh);
            ExtendedBidomainTissue<3> extended_tissue( &first_cell, &second_cell , &extra_factory);

            //set a value different from default for the conductivities of the second cell
            extended_tissue.SetIntracellularConductivitiesSecondCell(Create_c_vector(25.0,26.0,27.0));
            //this is normally done by the problem class, but in this test we do it manually
            extended_tissue.CreateIntracellularConductivityTensorSecondCell();
            extended_tissue.SetCacheReplication(cache_replication_saved); // Not the default to check it is archived...

            //shuffle default values to check if they get archived properly
            extended_tissue.SetAmFirstCell(11.0);
            extended_tissue.SetAmSecondCell(22.0);
            extended_tissue.SetAmGap(33.0);
            extended_tissue.SetCmFirstCell(44.0);
            extended_tissue.SetCmSecondCell(55.0);
            extended_tissue.SetGGap(66.0);

            //again, away from default value to check for archiving
            extended_tissue.SetUserSuppliedExtracellularStimulus(true);

            //set some heterogeneities in Ggap
            std::vector<boost::shared_ptr<AbstractChasteRegion<3> > > heterogeneity_areas;
            std::vector<double> Ggap_values;
            ChastePoint<3> cornerA(-1, -1, 0);
            ChastePoint<3> cornerB(0.001, 0.001, 0.001);
            boost::shared_ptr<ChasteCuboid<3> > p_cuboid_1(new ChasteCuboid<3>(cornerA, cornerB));
            heterogeneity_areas.push_back(p_cuboid_1);
            //within the first area
            Ggap_values.push_back(143.0);
            extended_tissue.SetGgapHeterogeneities(heterogeneity_areas, Ggap_values);
            extended_tissue.CreateGGapConductivities();

            // Some checks to make sure HeartConfig is being saved and loaded by this too.
            HeartConfig::Instance()->SetPrintingTimeStep(saved_printing_timestep);
            TS_ASSERT_DELTA(HeartConfig::Instance()->GetPrintingTimeStep(), saved_printing_timestep, 1e-9);

            intra_tensor_before_archiving = extended_tissue.rGetIntracellularConductivityTensor(0);
            intra_tensor_second_cell_before_archiving = extended_tissue.rGetIntracellularConductivityTensorSecondCell(0);
            extra_tensor_before_archiving = extended_tissue.rGetExtracellularConductivityTensor(0);

            // Save
            ArchiveOpener<boost::archive::text_oarchive, std::ofstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_oarchive* p_arch = arch_opener.GetCommonArchive();

            AbstractCardiacTissue<3>* const p_archive_bidomain_tissue = &extended_tissue;
            (*p_arch) << p_archive_bidomain_tissue;

            HeartConfig::Reset();
            TS_ASSERT_DELTA(HeartConfig::Instance()->GetPrintingTimeStep(), default_printing_timestep, 1e-9);
            TS_ASSERT_DIFFERS(saved_printing_timestep, default_printing_timestep);
        }
        //load
        {
            ArchiveOpener<boost::archive::text_iarchive, std::ifstream> arch_opener(archive_dir, archive_file);
            boost::archive::text_iarchive* p_arch = arch_opener.GetCommonArchive();

            AbstractCardiacTissue<3>* p_abstract_tissue;
            (*p_arch) >> p_abstract_tissue;
            assert(p_abstract_tissue!=NULL);

            //dynamic cast so we are able to test specific variables of ExtendedBidomainTissue
            ExtendedBidomainTissue<3>* p_extended_tissue = dynamic_cast<ExtendedBidomainTissue<3>*>(p_abstract_tissue);
            assert(p_extended_tissue != NULL);

            const c_matrix<double, 3, 3>& intra_tensor_after_archiving = p_extended_tissue->rGetIntracellularConductivityTensor(0);
//.........这里部分代码省略.........

    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));
//.........这里部分代码省略.........