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

    void TestWithHeterogeneousCellModels()
    {
        HeartConfig::Instance()->SetSimulationDuration(1.0); //ms
        HeartConfig::Instance()->SetUseStateVariableInterpolation(true);
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.005, 0.01, 1.0);

        TetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(0.01, 1.0);

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

        //// It's really very difficult to get hold of the correction assembler from here.
        //// The following, which requires 2 classes to be friends of this test compiles
        //// but fails asserts in the first line as bcc is not set up.
        //AbstractDynamicLinearPdeSolver<1,1,1>* p_solver = monodomain_problem.CreateSolver();
        //MonodomainSolver<1,1>* p_mono_solver = dynamic_cast<MonodomainSolver<1,1>*>(p_solver);
        //MonodomainCorrectionTermAssembler<1,1>* p_assembler = p_mono_solver->mpMonodomainCorrectionTermAssembler;
        //TS_ASSERT_EQUALS(p_assembler->mElementsCanDoSvi.size(), 10u);

        // Therefore, we just test that calling Solve() runs (without the checking that
        // cell models are identical, this fails).
        monodomain_problem.Solve();
    }

    void TestGetConductivityAndConductivityModifier() throw(Exception)
    {
        HeartConfig::Instance()->Reset();
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(1.0, 1.0); // [0,1] with h=1.0, ie 1 element mesh

        MyCardiacCellFactory<1> cell_factory;
        cell_factory.SetMesh(&mesh);

        BidomainTissue<1> bidomain_tissue( &cell_factory );

        double orig_extra_conductivity = bidomain_tissue.rGetExtracellularConductivityTensor(0)(0,0);
        double orig_intra_conductivity = bidomain_tissue.rGetIntracellularConductivityTensor(0)(0,0);

        TS_ASSERT_DELTA(orig_extra_conductivity, 7.0, 1e-9); // hard-coded using default
        TS_ASSERT_DELTA(orig_intra_conductivity, 1.75, 1e-9); // hard-coded using default

        SimpleConductivityModifier modifier;
        bidomain_tissue.SetConductivityModifier(&modifier);

        TS_ASSERT_DELTA(bidomain_tissue.rGetExtracellularConductivityTensor(0)(0,0), 2*orig_extra_conductivity, 1e-9);
        TS_ASSERT_DELTA(bidomain_tissue.rGetIntracellularConductivityTensor(0)(0,0), 2*orig_intra_conductivity, 1e-9);

        // The following asks for element 1 which doesn't exist
        TS_ASSERT_THROWS_THIS(bidomain_tissue.rGetExtracellularConductivityTensor(1)(0,0),
                              "Conductivity tensor requested for element with global_index=1, but there are only 1 elements in the mesh.")
    }

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

    void TestSolveCellSystemsInclUpdateVoltage() throw(Exception)
    {
        HeartConfig::Instance()->Reset();
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(1.0, 1.0); // [0,1] with h=1.0, ie 2 node mesh

        MyCardiacCellFactory cell_factory;
        cell_factory.SetMesh(&mesh);

        MonodomainTissue<1> monodomain_tissue( &cell_factory );

        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))
        {
            TS_ASSERT_DELTA(monodomain_tissue.GetCardiacCell(0)->GetVoltage(), -81.4354, 1e-3);
        }

        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(1))
        {
            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))
        {
            // 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))
        {
            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);
        }

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

    void TestCoverage3d()
    {
        HeartConfig::Instance()->SetSimulationDuration(0.1); //ms
        HeartConfig::Instance()->SetUseStateVariableInterpolation(true);
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.005, 0.01, 0.1);

        TetrahedralMesh<3,3> mesh;
        mesh.ConstructRegularSlabMesh(0.02, 0.02, 0.02, 0.02);

        ZeroStimulusCellFactory<CellLuoRudy1991FromCellML,3> cell_factory;
        MonodomainProblem<3> monodomain_problem( &cell_factory );
        monodomain_problem.SetMesh(&mesh);
        monodomain_problem.Initialise();
        monodomain_problem.Solve();
    }

    void RunBenchMark(double h, double dt, double endTime, bool useSvi)
    {
        TetrahedralMesh<3,3> mesh;

        mesh.ConstructRegularSlabMesh(h, 2.0, 0.7, 0.3);

        std::stringstream output_dir;
        output_dir << "Benchmark" << "_h" << h << "_dt" << dt;

        HeartConfig::Instance()->SetOutputDirectory(output_dir.str());
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");
        HeartConfig::Instance()->SetSimulationDuration(endTime); //ms

        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.005, dt, 0.1);
        HeartConfig::Instance()->SetSurfaceAreaToVolumeRatio(1400); // 1400 1/cm
        HeartConfig::Instance()->SetCapacitance(1); // 1uF/cm^2

        HeartConfig::Instance()->SetVisualizeWithMeshalyzer(false);

        // The Chaste results for the benchmark paper use STATE-VARIABLE INTERPOLATION switched on
        // (see comments above)
        HeartConfig::Instance()->SetUseStateVariableInterpolation(useSvi);

        // Regarding the second paper described above, to run the simulations with ICI, comment out the
        // above line. To run the simulation with operator splitting, or with (full) mass-lumping,
        // comment out the above SVI line and uncomment one of the below. (Note: half-lumping is not
        // available).
        //HeartConfig::Instance()->SetUseMassLumping(true); // what is described as full-lumping in this paper
        //HeartConfig::Instance()->SetUseReactionDiffusionOperatorSplitting(true);

        double long_conductance = 0.17 * 0.62/(0.17+0.62) * 10;    // harmonic mean of 0.17, 0.62 S/m converted to mS/cm
        double trans_conductance = 0.019 * 0.24/(0.019+0.24) * 10; // harmonic mean of 0.019,0.24 S/m converted to mS/cm

        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(long_conductance, trans_conductance, trans_conductance));

        BenchmarkCellFactory cell_factory;

        MonodomainProblem<3> problem( &cell_factory );
        problem.SetMesh(&mesh);

        problem.Initialise();
        problem.SetWriteInfo();
        problem.Solve();
    }

    void TestMonodomainTissueGetCardiacCell() throw(Exception)
    {
        HeartConfig::Instance()->Reset();
        TetrahedralMesh<1,1> mesh;
        mesh.ConstructRegularSlabMesh(1.0, 1.0); // [0,1] with h=1.0, ie 2 node mesh

        MyCardiacCellFactory cell_factory;
        cell_factory.SetMesh(&mesh);

        MonodomainTissue<1> monodomain_tissue( &cell_factory );

        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(0))
        {
            AbstractCardiacCellInterface* cell = monodomain_tissue.GetCardiacCell(0);
            TS_ASSERT_DELTA(cell->GetStimulus(0.001),-80,1e-10);
        }

        if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(1))
        {
            AbstractCardiacCellInterface* cell = monodomain_tissue.GetCardiacCell(1);
            TS_ASSERT_DELTA(cell->GetStimulus(0.001),0,1e-10);
        }
    }

    void TestMonodomain3d() throw(Exception)
    {
        /* HOW_TO_TAG Cardiac/Problem definition
         * Generate a slab (cuboid) mesh rather than read a mesh in, and pass it to solver
         */

        /* We will auto-generate a mesh this time, and pass it in, rather than
         * provide a mesh file name. This is how to generate a cuboid mesh with
         * a given spatial stepsize h */
        TetrahedralMesh<3,3> mesh;
        double h=0.02;
        mesh.ConstructRegularSlabMesh(h, 0.8 /*length*/, 0.3 /*width*/, 0.3 /*depth*/);
        /* (In 2D the call is identical, but without the depth parameter).
         *
         * EMPTYLINE
         *
         * Set the simulation duration, etc, and create an instance of the cell factory.
         * One thing that should be noted for monodomain problems, the ''intracellular
         * conductivity'' is used as the monodomain effective conductivity (not a
         * harmonic mean of intra and extracellular conductivities). So if you want to
         * alter the monodomain conductivity call
         * `HeartConfig::Instance()->SetIntracellularConductivities`
         */
        HeartConfig::Instance()->SetSimulationDuration(5); //ms
        HeartConfig::Instance()->SetOutputDirectory("Monodomain3dExample");
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");
        HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.005, 0.01, 0.1);

        BenchmarkCellFactory cell_factory;

        /* Now we declare the problem class, `MonodomainProblem<3>` instead of `BidomainProblem<2>`.
         * The interface for both is the same.
         */
        MonodomainProblem<3> monodomain_problem( &cell_factory );

        /* If a mesh-file-name hasn't been set using `HeartConfig`, we have to pass in
         * a mesh using the `SetMesh` method (must be called before `Initialise`). */
        monodomain_problem.SetMesh(&mesh);

        /* By default data for all nodes is output, but for big simulations, sometimes this
         * might not be required, and the action potential only at certain nodes required.
         * The following code shows how to output the results at the first, middle and last
         * nodes, for example. (The output is written to the HDF5 file; regular visualisation output
         * will be turned off. HDF5 files can be read using Matlab). We are not using this in this
         * simulation however (hence the boolean being set to false).
         */
        bool partial_output = false;
        if(partial_output)
        {
            std::vector<unsigned> nodes_to_be_output;
            nodes_to_be_output.push_back(0);
            nodes_to_be_output.push_back((unsigned)round( (mesh.GetNumNodes()-1)/2 ));
            nodes_to_be_output.push_back(mesh.GetNumNodes()-1);
            monodomain_problem.SetOutputNodes(nodes_to_be_output);
        }

        /* `SetWriteInfo` is a useful method that means that the min/max voltage is
         * printed as the simulation runs (useful for verifying that cells are stimulated
         * and the wave propagating, for example) (although note scons does buffer output
         * before printing to screen) */
        monodomain_problem.SetWriteInfo();

        /* Finally, call `Initialise` and `Solve` as before */
        monodomain_problem.Initialise();
        monodomain_problem.Solve();

        /* This part is just to check nothing has accidentally been changed in this example */
        ReplicatableVector voltage(monodomain_problem.GetSolution());
        TS_ASSERT_DELTA(voltage[0], 34.9032, 1e-2);
    }

    /*
     * == Test 2: Solving a linear parabolic PDE ==
     *
     * Now we solve a parabolic PDE. We choose a simple problem so that the code changes
     * needed from the elliptic case are clearer. We will solve
     * du/dt = div(grad u) + u, in 3d, with boundary conditions u=1 on the boundary, and initial
     * conditions u=1.
     *
     */
    void TestSolvingParabolicPde() throw(Exception)
    {
        /* Create a 10 by 10 by 10 mesh in 3D, this time using the {{{ConstructRegularSlabMesh}}} method
         * on the mesh. The first parameter is the cartesian space-step and the other three parameters are the width, height and depth of the mesh.*/
        TetrahedralMesh<3,3> mesh;
        mesh.ConstructRegularSlabMesh(0.1, 1.0, 1.0, 1.0);

        /* Our PDE object should be a class that is derived from the {{{AbstractLinearParabolicPde}}}.
         * We could write it ourselves as in the previous test, but since the PDE we want to solve is
         * so simple, it has already been defined (look it up! - it is located in pde/test/pdes).
         */
        HeatEquationWithSourceTerm<3> pde;

        /* Create a new boundary conditions container and specify u=1.0 on the boundary. */
        BoundaryConditionsContainer<3,3,1> bcc;
        bcc.DefineConstantDirichletOnMeshBoundary(&mesh, 1.0);

        /* Create an instance of the solver, passing in the mesh, pde and boundary conditions.
         */
        SimpleLinearParabolicSolver<3,3> solver(&mesh,&pde,&bcc);

        /* For parabolic problems, initial conditions are also needed. The solver will expect
         * a PETSc vector, where the i-th entry is the initial solution at node i, to be passed
         * in. To create this PETSc {{{Vec}}}, we will use a helper function in the {{{PetscTools}}}
         * class to create a {{{Vec}}} of size num_nodes, with each entry set to 1.0. Then we
         * set the initial condition on the solver. */
        Vec initial_condition = PetscTools::CreateAndSetVec(mesh.GetNumNodes(), 1.0);
        solver.SetInitialCondition(initial_condition);

        /* Next define the start time, end time, and timestep, and set them. */
        double t_start = 0;
        double t_end = 1;
        double dt = 0.01;
        solver.SetTimes(t_start, t_end);
        solver.SetTimeStep(dt);


        /* HOW_TO_TAG PDE
         * Output results to file for time-dependent PDE solvers
         */


        /* When we call Solve() below we will just get the solution at the final time. If we want
         * to have intermediate solutions written to file, we do the following. We start by
         * specifying an output directory and filename prefix for our results file:
         */
        solver.SetOutputDirectoryAndPrefix("ParabolicSolverTutorial","results");
        /* When an output directory has been specified, the solver writes output in HDF5 format. To
         * convert this to another output format, we call the relevant method. Here, we convert
         * the output to plain text files (VTK or cmgui formats are also possible). We also say how
         * often to write the data, telling the solver to output results to file every tenth timestep.
         * The solver will create one file for each variable (in this case there is only one variable)
         * and for each time, so for example, the file
         * `results_Variable_0_10` is the results for u, over all nodes, at the 11th printed time.
         * Have a look in the output directory after running the test. (For comments on visualising the data in
         * matlab or octave, see the end of the tutorial UserTutorials/WritingPdeSolvers.)
         */
        solver.SetOutputToTxt(true);
        solver.SetPrintingTimestepMultiple(10);

        /* Now we can solve the problem. The {{{Vec}}} that is returned can be passed into a
         * {{{ReplicatableVector}}} as before.
         */
        Vec solution = solver.Solve();
        ReplicatableVector solution_repl(solution);

        /* Let's also solve the equivalent static PDE, i.e. set du/dt=0, so 0=div(gradu) + u. This
         * is easy, as the PDE class has already been defined. */
        SimplePoissonEquation<3,3> static_pde;
        SimpleLinearEllipticSolver<3,3> static_solver(&mesh, &static_pde, &bcc);
        Vec static_solution = static_solver.Solve();
        ReplicatableVector static_solution_repl(static_solution);

        /* We can now compare the solution of the parabolic PDE at t=1 with the static solution,
         * to see if the static equilibrium solution was reached in the former. (Ideally we should
         * compute some relative error, but we just compute an absolute error for simplicity.) */
        for (unsigned i=0; i<static_solution_repl.GetSize(); i++)
        {
            TS_ASSERT_DELTA( solution_repl[i], static_solution_repl[i], 1e-3);
        }

        /* All PETSc vectors should be destroyed when they are no longer needed. */
        PetscTools::Destroy(initial_condition);
        PetscTools::Destroy(solution);
        PetscTools::Destroy(static_solution);
    }

    // this method loads the output file from the previous method and computes the activation
    // time (defined as the time V becomes positive) for each node.
    void ConvertToActivationMap(double h, double dt, bool useSvi)
    {
        //TetrahedralMesh<3,3> mesh1;
        //double h1=0.01;    // 0.01, 0.02, 0.05
        //mesh1.ConstructRegularSlabMesh(h1, 2.0, 0.7, 0.3);
        //MeshalyzerMeshWriter<3,3> writer("Mesh0.01", "mesh01");
        //writer.WriteFilesUsingMesh(mesh1);

        TetrahedralMesh<3,3> mesh;
        double printing_dt=0.1;
        mesh.ConstructRegularSlabMesh(h, 2.0, 0.7, 0.3);

        std::stringstream input_dir;
        input_dir << "Benchmark" << "_h" << h << "_dt" << dt;
        Hdf5DataReader reader(input_dir.str(),"results");

        unsigned num_timesteps = reader.GetUnlimitedDimensionValues().size();
        DistributedVectorFactory factory(mesh.GetNumNodes());
        Vec voltage = factory.CreateVec();


        std::vector<double> activation_times(mesh.GetNumNodes(), -1.0);
        std::vector<double> last_negative_voltage(mesh.GetNumNodes(), 1.0);
        for(unsigned timestep=0; timestep<num_timesteps; timestep++)
        {
            reader.GetVariableOverNodes(voltage, "V", timestep);
            ReplicatableVector voltage_repl(voltage);

            for(unsigned i=0; i<mesh.GetNumNodes(); i++)
            {
                double V = voltage_repl[i];
                if(V > 0 && activation_times[i] < 0.0)
                {
                    double old = last_negative_voltage[i];
                    assert(old < 0);
                    activation_times[i] = (timestep-V/(V-old))*printing_dt;
                }
                else if (V<=0)
                {
                    last_negative_voltage[i]=V;
                }
            }
        }

        OutputFileHandler handler("ActivationMaps", false);
        if (PetscTools::AmMaster() == false)
        {
            return;
        }
        //Only master proceeds to write


        c_vector<double, 3> top_corner;
        top_corner[0] = 2.0;
        top_corner[1] = 0.7;
        top_corner[2] = 0.3;
        c_vector<double, 3> unit_diagonal = top_corner/norm_2(top_corner);

        std::stringstream output_file1;
        output_file1 << "diagonal" << "_h" << h << "_dt" << dt;
        if (useSvi)
        {
            output_file1 << "_svi.dat";
        }
        else
        {
            output_file1 << "_ici.dat";
        }
        out_stream p_diag_file = handler.OpenOutputFile(output_file1.str());

        for(unsigned i=0; i<mesh.GetNumNodes(); i++)
        {
            c_vector<double, 3> position =  mesh.GetNode(i)->rGetLocation();
            c_vector<double, 3>  projected_diagonal = unit_diagonal*inner_prod(unit_diagonal, position);
            double off_diagonal = norm_2(position - projected_diagonal);

            if (off_diagonal < h/3)
            {
                double distance = norm_2(position);
                (*p_diag_file) << distance<<"\t"<< activation_times[i]<<"\t"<<off_diagonal<<"\n";
                if( fabs(position[0]-2.0) < 1e-8)
                {
                    std::cout << "h, dt = " << h << ", " << dt << "\n\t";
                    std::cout << "activation_times[" << i << "] = " << activation_times[i] << "\n";
                }
            }
        }
        p_diag_file->close();

        std::stringstream output_file;
        output_file << "activation" << "_h" << h << "_dt" << dt << ".dat";
        out_stream p_file = handler.OpenOutputFile(output_file.str());

        for(unsigned i=0; i<activation_times.size(); i++)
        {
            *p_file << activation_times[i] << "\n";
        }
        p_file->close();
//.........这里部分代码省略.........

    // Test the functionality for outputting the values of requested cell state variables
    void TestExtendedBidomainProblemPrintsMultipleVariables() throw (Exception)
    {
        // Get the singleton in a clean state
        HeartConfig::Instance()->Reset();

        // Set configuration file
        std::string dir = "ExtBidoMultiVars";
        std::string filename = "extended";
        HeartConfig::Instance()->SetOutputDirectory(dir);
        HeartConfig::Instance()->SetOutputFilenamePrefix(filename);

        HeartConfig::Instance()->SetSimulationDuration(0.1);
        // HeartConfig::Instance()->SetKSPSolver("gmres");
        HeartConfig::Instance()->SetUseAbsoluteTolerance(2e-4);
        HeartConfig::Instance()->SetKSPPreconditioner("jacobi");

        /** Check that also the converters handle multiple variables**/
        HeartConfig::Instance()->SetVisualizeWithCmgui(true);
        HeartConfig::Instance()->SetVisualizeWithMeshalyzer(true);

        TetrahedralMesh<1,1> mesh;
        unsigned number_of_elements = 100;
        double length = 10.0;
        mesh.ConstructRegularSlabMesh(length/number_of_elements, length);

        // Override the variables we are interested in writing.
        std::vector<std::string> output_variables;
        output_variables.push_back("calcium_dynamics__Ca_NSR");
        output_variables.push_back("ionic_concentrations__Nai");
        output_variables.push_back("fast_sodium_current_j_gate__j");
        output_variables.push_back("ionic_concentrations__Ki");

        HeartConfig::Instance()->SetOutputVariables( output_variables );

        // Set up problem
        PlaneStimulusCellFactory<CellFaberRudy2000FromCellML, 1> cell_factory_1(-60, 0.5);
        PlaneStimulusCellFactory<CellLuoRudy1991FromCellML, 1> cell_factory_2(0.0,0.5);
        ExtendedBidomainProblem<1> ext_problem( &cell_factory_1, &cell_factory_2 );
        ext_problem.SetMesh(&mesh);

        // Solve
        ext_problem.Initialise();
        ext_problem.Solve();

        // Get a reference to a reader object for the simulation results
        Hdf5DataReader data_reader1 = ext_problem.GetDataReader();
        std::vector<double> times = data_reader1.GetUnlimitedDimensionValues();

        // Check there is information about 11 timesteps (0, 0.01, 0.02, ...)
        unsigned num_steps = 11u;
        TS_ASSERT_EQUALS( times.size(), num_steps);
        TS_ASSERT_DELTA( times[0], 0.0, 1e-12);
        TS_ASSERT_DELTA( times[1], 0.01, 1e-12);
        TS_ASSERT_DELTA( times[2], 0.02, 1e-12);
        TS_ASSERT_DELTA( times[3], 0.03, 1e-12);

        // There should be 11 values per variable and node.
        std::vector<double> node_5_v = data_reader1.GetVariableOverTime("V", 5);
        TS_ASSERT_EQUALS( node_5_v.size(), num_steps);

        std::vector<double> node_5_v_2 = data_reader1.GetVariableOverTime("V_2", 5);
        TS_ASSERT_EQUALS( node_5_v_2.size(), num_steps);

        std::vector<double> node_5_phi = data_reader1.GetVariableOverTime("Phi_e", 5);
        TS_ASSERT_EQUALS( node_5_phi.size(), num_steps);

        for (unsigned i=0; i<output_variables.size(); i++)
        {
            unsigned global_index = 2+i*2;
            std::vector<double> values = data_reader1.GetVariableOverTime(output_variables[i], global_index);
            TS_ASSERT_EQUALS( values.size(), num_steps);

            // Check the last values match the cells' state
            if (ext_problem.rGetMesh().GetDistributedVectorFactory()->IsGlobalIndexLocal(global_index))
            {
                AbstractCardiacCellInterface* p_cell = ext_problem.GetTissue()->GetCardiacCell(global_index);
                TS_ASSERT_DELTA(values.back(), p_cell->GetAnyVariable(output_variables[i],0), 1e-12);
            }

            //check the extra files for extra variables are there (the content is tested in the converter's tests)
            FileFinder file(dir + "/output/"+ filename +"_"+ output_variables[i] + ".dat", RelativeTo::ChasteTestOutput);
            TS_ASSERT(file.Exists());
        }
    }

    /**
     * This test is aimed at comparing the extended bidomain implementation in Chaste with
     * the original Finite Difference code developed by Martin Buist.
     *
     * All the parameters are chosen to replicate the same conditions as in his code.
     */
    void TestExtendedProblemVsMartincCode() throw (Exception)
    {
        SetupParameters();

        TetrahedralMesh<1,1> mesh;
        unsigned number_of_elements = 100;//this is nGrid in Martin's code
        double length = 10.0;//100mm as in Martin's code
        mesh.ConstructRegularSlabMesh(length/number_of_elements, length);
        TS_ASSERT_EQUALS(mesh.GetNumAllNodes(), number_of_elements + 1);

        double Am_icc = 1000.0;
        double Am_smc = 1000.0;
        double Am_gap = 1.0;
        double Cm_icc = 1.0;
        double Cm_smc = 1.0;
        double G_gap = 20.0;//mS/cm^2

        HeartConfig::Instance()->SetSimulationDuration(1000.0);  //ms.

        ICC_Cell_factory icc_factory;
        SMC_Cell_factory smc_factory;

        std::string dir = "ICCandSMC";
        std::string filename = "extended1d";
        HeartConfig::Instance()->SetOutputDirectory(dir);
        HeartConfig::Instance()->SetOutputFilenamePrefix(filename);

        ExtendedBidomainProblem<1> extended_problem( &icc_factory , &smc_factory);
        extended_problem.SetMesh(&mesh);

        extended_problem.SetExtendedBidomainParameters(Am_icc,Am_smc, Am_gap, Cm_icc, Cm_smc, G_gap);
        extended_problem.SetIntracellularConductivitiesForSecondCell(Create_c_vector(1.0));

        std::vector<unsigned> outputnodes;
        outputnodes.push_back(50u);
        HeartConfig::Instance()->SetRequestedNodalTimeTraces(outputnodes);

        extended_problem.Initialise();
        extended_problem.Solve();
        HeartEventHandler::Headings();
        HeartEventHandler::Report();

        /**
         * Compare with valid data.
         * As Martin's code is an FD code, results will never match exactly.
         * The comparison below is done against a 'valid' h5 file.
         *
         * The h5 file (1DValid.h5) is a Chaste (old phi_i formulation) file with is valid because, when extrapolating results from it, they look very similar
         * (except for a few points at the end of the upstroke) to the results taken
         * directly from Martin's code.
         * A plot of Chaste results versus Martin's result (at node 50) is stored
         * in the file 1DChasteVsMartin.eps for reference.
         *
         * A second plot comparing the old formulation (with phi_i) to the new formulation with V_m is contained in
         *.1DChasteNewFormulation.png
         *
         */
         TS_ASSERT( CompareFilesViaHdf5DataReader("heart/test/data/extendedbidomain", "1DValid", false,
                                                dir, filename, true,
                                                0.2));
         /*
          * Here we compare the new formulation (V_m1, V_m2, phi_e)
          *  with the previous formulation (phi_i1, phi_i2, phi_e) running with GMRES and an absolute KSP tolerance of 1e-8.
          */
         TS_ASSERT( CompareFilesViaHdf5DataReader("heart/test/data/extendedbidomain", "extended1d_previous_chaste_formulation_abs_tol_1e-8", false,
                                                dir, filename, true,
                                                1e-2));
    }

    // like TestOperatorSplittingMonodomainSolver but much finer mesh and smaller
    // dt so can check for proper convergence
    void TestConvergenceOnFineMesh() throw(Exception)
    {
        ReplicatableVector final_voltage_normal;
        ReplicatableVector final_voltage_operator_splitting;

        HeartConfig::Instance()->SetSimulationDuration(4.0); //ms
        HeartConfig::Instance()->SetOutputFilenamePrefix("results");

        double h = 0.001; // very fine

        // Normal
        {
            HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.001, 0.001, 0.1); // very small timesteps

            TetrahedralMesh<1,1> mesh;
            mesh.ConstructRegularSlabMesh(h, 1.0);
            HeartConfig::Instance()->SetOutputDirectory("MonodomainCompareWithOperatorSplitting_normal");
            BlockCellFactory<1> cell_factory;

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

            final_voltage_normal.ReplicatePetscVector(monodomain_problem.GetSolution());
        }

        // Operator splitting
        {
            HeartConfig::Instance()->SetOdePdeAndPrintingTimeSteps(0.001, 0.002, 0.1); // very small timesteps - use pde-dt = 2 ode-dt
                                                                                       // as effective_ode_dt = min{pde_dt/2, ode_dt} in
                                                                                       // our operator splitting implementation

            TetrahedralMesh<1,1> mesh;
            mesh.ConstructRegularSlabMesh(h, 1.0);
            HeartConfig::Instance()->SetOutputDirectory("MonodomainCompareWithOperatorSplitting_splitting");
            BlockCellFactory<1> cell_factory;

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

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

            final_voltage_operator_splitting.ReplicatePetscVector(monodomain_problem.GetSolution());
        }

        bool some_node_depolarised = false;
        assert(final_voltage_normal.GetSize()==final_voltage_operator_splitting.GetSize());
        for(unsigned j=0; j<final_voltage_normal.GetSize(); j++)
        {
            double tol=4.7;

            TS_ASSERT_DELTA(final_voltage_normal[j], final_voltage_operator_splitting[j], tol);

            if(final_voltage_normal[j]>-80)
            {
                // shouldn't be exactly equal, as long as away from resting potential
                TS_ASSERT_DIFFERS(final_voltage_normal[j], final_voltage_operator_splitting[j]);
            }

            if(final_voltage_normal[j]>0.0)
            {
                some_node_depolarised = true;
            }
        }
        UNUSED_OPT(some_node_depolarised);
        assert(some_node_depolarised);
    }

    void TestConductionVelocityInCrossFibreDirection2d() throw(Exception)
    {
        ReplicatableVector final_solution_ici;
        ReplicatableVector final_solution_svi;

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

        // much lower conductivity in cross-fibre direction - ICI will struggle
        HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(1.75, 0.17));
        HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(7.0, 0.7));

        // ICI - nodal current interpolation - the default
        {
            TetrahedralMesh<2,2> mesh;
            mesh.ConstructRegularSlabMesh(0.02 /*h*/, 0.5, 0.3);

            HeartConfig::Instance()->SetOutputDirectory("BidomainIci2d");
            HeartConfig::Instance()->SetOutputFilenamePrefix("results");

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

            BlockCellFactory<2> cell_factory;
            BidomainProblem<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
        {
            TetrahedralMesh<2,2> mesh;
            mesh.ConstructRegularSlabMesh(0.02 /*h*/, 0.5, 0.3);

            HeartConfig::Instance()->SetOutputDirectory("BidomainSvi2d");
            HeartConfig::Instance()->SetOutputFilenamePrefix("results");

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

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

            bidomain_problem.Solve();

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

        // See comments in equivalent part of test/monodomain/TestMonodomainWithSvi.hpp
        // For bidomain with h=0.02, SVI CV is slower in both fibre and cross-fibre
        // directions

        // node 20 (for h=0.02) is on the x-axis (fibre direction)
        TS_ASSERT_DELTA(final_solution_ici[20*2], -64.1105, 1e-3); // Node 20 phi_i
        TS_ASSERT_DELTA(final_solution_svi[20*2], -78.0936, 1e-3); // Node 20 phi_i
        // node 234 (for h=0.02) is on the y-axis (cross-fibre direction)
        TS_ASSERT_DELTA(final_solution_ici[234*2], -57.7239, 1e-3); // Node 234 phi_i
        TS_ASSERT_DELTA(final_solution_svi[234*2], 38.9004, 1e-3); // Node 234 phi_i
    }

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