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

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

//.........这里部分代码省略.........
            {
                for(unsigned j=0; j<3; j++)
                {
                    TS_ASSERT_DELTA(intra_tensor_before_archiving(i,j), intra_tensor_after_archiving(i,j), 1e-9);
                    TS_ASSERT_DELTA(intra_tensor_second_cell_before_archiving(i,j), intra_tensor_second_cell_after_archiving(i,j), 1e-9);
                    TS_ASSERT_DELTA(extra_tensor_before_archiving(i,j), extra_tensor_after_archiving(i,j), 1e-9);
                }
            }

            //check that the member variable mIntracellularConductivitiesSecondCell was archived properly
            TS_ASSERT_EQUALS(p_extended_tissue->GetIntracellularConductivitiesSecondCell()(0),25.0);
            TS_ASSERT_EQUALS(p_extended_tissue->GetIntracellularConductivitiesSecondCell()(1),26.0);
            TS_ASSERT_EQUALS(p_extended_tissue->GetIntracellularConductivitiesSecondCell()(2),27.0);

            //check that we get the same values from the archive which are different from the default
            TS_ASSERT_EQUALS(p_extended_tissue->GetAmFirstCell(), 11.0);
            TS_ASSERT_EQUALS(p_extended_tissue->GetAmSecondCell(), 22.0);
            TS_ASSERT_EQUALS(p_extended_tissue->GetAmGap(), 33.0);
            TS_ASSERT_EQUALS(p_extended_tissue->GetCmFirstCell(), 44.0);
            TS_ASSERT_EQUALS(p_extended_tissue->GetCmSecondCell(), 55.0);
            TS_ASSERT_EQUALS(p_extended_tissue->GetGGap(), 66.0);

            // We shouldn't need to re-build the mesh, but we use it to check that the new tissue has the same mesh
            // Also, when testing in parallel, we use it to get the vector factory to loop over the nodes we own.
            // this is because  p_extended_tissue->pGetMesh()->GetDistributedVectorFactory() doesn't compile (discards qualifier stuff caused by use of const).
            TrianglesMeshReader<3,3> mesh_reader("mesh/test/data/cube_136_elements");
            DistributedTetrahedralMesh<3,3> mesh;
            mesh.ConstructFromMeshReader(mesh_reader);

            TS_ASSERT_EQUALS(mesh.GetNumNodes(), p_extended_tissue->pGetMesh()->GetNumNodes());//note: this is allowed because GetNumNodes has const in the signature

            //check archiving of stimulus for first cell at some random times (it is unstimulated everywhere at all times)
            for (unsigned i = 0; i < mesh.GetNumNodes(); i++)
            {
                if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(i))
                {
                    TS_ASSERT_EQUALS(p_extended_tissue->GetCardiacCell(i)->GetIntracellularStimulus(0.0), 0.0);
                    TS_ASSERT_EQUALS(p_extended_tissue->GetCardiacCell(i)->GetIntracellularStimulus(0.1), 0.0);
                    TS_ASSERT_EQUALS(p_extended_tissue->GetCardiacCell(i)->GetIntracellularStimulus(2.5), 0.0);
                    TS_ASSERT_EQUALS(p_extended_tissue->GetCardiacCell(i)->GetIntracellularStimulus(28.9), 0.0);
                }
            }

            //for second cell and other stuff, we probe nodes 0 and 1.
            unsigned node_0 = 0u;
            unsigned node_1 = 1u;
            //If the test is run in parallel, we need to work out the new indices
            const std::vector<unsigned>& r_permutation = mesh.rGetNodePermutation();
            if (!r_permutation.empty())
            {
                node_0 = r_permutation[0u];
                node_1 = r_permutation[1u];
            }

            //second cell is stimulated in the corner (node 0) from time 0 to 1. check it gets all this after loading
            if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(node_0))
            {
                TS_ASSERT_EQUALS(p_extended_tissue->GetCardiacSecondCell(node_0)->GetIntracellularStimulus(0.5), -105.0*1400);
                TS_ASSERT_EQUALS(p_extended_tissue->GetCardiacSecondCell(node_0)->GetIntracellularStimulus(2.5), 0.0);

                //find local index of (the new) node_0, it should be in the heterogeneity region
                unsigned ownership_range_low = mesh.GetDistributedVectorFactory()->GetLow();
                unsigned  local_index = node_0 - ownership_range_low;
                //std::cout<<local_index<<std::endl;
                TS_ASSERT_EQUALS(p_extended_tissue->rGetGapsDistributed()[local_index],143.0);//g_gap value inside heterogeneity region
            }

            //node 0 has extracellular stimulus (1 ms from 0.1)
            if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(node_0))
            {
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_0)->GetStimulus(0.0), 0);
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_0)->GetStimulus(0.5), -428000);
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_0)->GetStimulus(1.0), -428000);
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_0)->GetStimulus(1.15), 0);
            }

            //node 1 doesn't
            if (mesh.GetDistributedVectorFactory()->IsGlobalIndexLocal(node_1))
            {
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_1)->GetStimulus(0.0), 0);
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_1)->GetStimulus(0.5), 0);
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_1)->GetStimulus(1.0), 0);
                TS_ASSERT_EQUALS(p_extended_tissue->GetExtracellularStimulus(node_1)->GetStimulus(1.15), 0);

                //find local index of (the new) node_1, it should NOT be in the heterogeneity region
                unsigned ownership_range_low = mesh.GetDistributedVectorFactory()->GetLow();
                unsigned  local_index = node_1 - ownership_range_low;
                TS_ASSERT_EQUALS(p_extended_tissue->rGetGapsDistributed()[local_index],66.0);//standard g_gap value, outside heterogeneity region
            }

            //check the archiving of the flag (it would be false by default, but we set it to true before archiving)
            TS_ASSERT_EQUALS(p_extended_tissue->HasTheUserSuppliedExtracellularStimulus(),true);

            TS_ASSERT_EQUALS(cache_replication_saved, p_extended_tissue->GetDoCacheReplication());
            TS_ASSERT_DELTA(HeartConfig::Instance()->GetPrintingTimeStep(), saved_printing_timestep, 1e-9);
            TS_ASSERT_DIFFERS(saved_printing_timestep, default_printing_timestep); // Test we are testing something in case default changes

            delete p_extended_tissue;
        }
    }

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