Another Bidomain Simulation
This tutorial is automatically generated from TestAnotherBidomainSimulationTutorial.hpp at revision 3c544f98da9c. Note that the code is given in full at the bottom of the page.
Another example showing how to run a bidomain simulation
In this tutorial we run another bidomain simulation, showing:
- i) an example using one of the source cell factories,
- ii) how to define and use fibre directions, and
- iii) mentioning how to write other output file formats.
The first thing to do is to include the headers as before.
#include <cxxtest/TestSuite.h>
#include "BidomainProblem.hpp"
#include "LuoRudy1991.hpp"
#include "PetscSetupAndFinalize.hpp"
Rather than write our own cell factory this time, we will use
the PlaneStimulusCellFactory
.
#include "PlaneStimulusCellFactory.hpp"
Now we define the test class, which must inherit from CxxTest::TestSuite
as usual, and the (public) test method
class TestAnotherBidomainTutorial : public CxxTest::TestSuite
{
public:
void TestFibreSimulation()
{
It is not the case here, but if there were other tests in the file that
had already been run and might have changed parameters in HeartConfig
, we
would need to call Reset
HeartConfig::Instance()->Reset();
Next, we have to create a cell factory of the type we defined above. The plane
stimulus cell factory sets up cells of the given type with a non-zero stimulus
for cells on the $x=0$ boundary. The 2 below is the dimension, and the -2000000
is the stimulus magnitude (the default stimulus duration is used).
PlaneStimulusCellFactory<CellLuoRudy1991FromCellML,2> cell_factory(-2000000);
Define an end time, output directory and prefix as before
HeartConfig::Instance()->SetSimulationDuration(5.0); //ms
HeartConfig::Instance()->SetOutputDirectory("BidomainFibresTutorial");
HeartConfig::Instance()->SetOutputFilenamePrefix("results");
Define a mesh to be read, saying that we also want to read fibres. The extra part can either be
cp::media_type::Orthotropic
, in which case 2D_0_to_1mm_800_elements.ortho
will also be read;
or cp::media_type::Axisymmetric
, in which case 2D_0_to_1mm_800_elements.axi
will also be read.
See the file formats documentation for full descriptions of these formats, but basically .axi
files provide the fibre direction for each element in the mesh, and .ortho files provide the fibre,
sheet (and normal in 3D) directions for each element in the mesh.
HeartConfig::Instance()->SetMeshFileName("mesh/test/data/2D_0_to_1mm_800_elements", cp::media_type::Orthotropic);
The fibre file provided here defines (non-physiological) ‘kinked’ fibres which are
in the x-direction for $x<0.05$, and then the diagonal $(1,1)$ direction for $x \geqslant 0.05$.
It was generated by looping over the mesh that is to be used, as shown here.
The output is the .ortho
file minus the header.
// TetrahedralMesh<2,2> mesh;
// TrianglesMeshReader<2,2> reader("mesh/test/data/2D_0_to_1mm_800_elements");
// mesh.ConstructFromMeshReader(reader);
// for (unsigned i=0; i<mesh.GetNumElements(); i++)
// {
// double x = mesh.GetElement(i)->CalculateCentroid()[0];
// if (x<0.05)
// {
// std::cout << "1 0 0 1\n";
// }
// else
// {
// std::cout << "0.707106781 0.707106781 -0.707106781 0.707106781\n";
// }
// }
// return;
This is of course not enough - by default isotropic conductivities are used so the variable fibre directions won’t make any difference - so we have to alter the intracellular and extracellular conductivities to be anisotropic. The first value in the following is the conductivity in the fibre direction, the second the conductivity in the sheet direction (and the third in 3d would be that in the normal direction). Note, we have scaled all the conductivities from the physiological values as the mesh is a little too coarse to be able to handle the smallest conductivities (try running with scale = 1 to see the error message).
double scale = 2;
HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(1.75*scale, 0.19*scale));
HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(7.0*scale, 2.4*scale));
The output will be written to $CHASTE_TEST_OUTPUT/BidomainFibresTutorial
in HDF5 format. It can be converted to visualisable formats (Meshalyzer, Cmgui or VTK) at the end
of the simulation using methods in HeartConfig
, e.g.
//HeartConfig::Instance()->SetVisualizeWithMeshalyzer(true);
HeartConfig::Instance()->SetVisualizeWithCmgui(true);
The other option is to write in VTK format (which needs VTK installed), following which the results can be loaded in the visualiser Paraview
//HeartConfig::Instance()->SetVisualizeWithVtk(true);
If the mesh is a DistributedTetrahedralMesh then we can use parallel VTK files (.pvtu)
//HeartConfig::Instance()->SetVisualizeWithParallelVtk(true);
Now we create a problem class, initialise and solve
BidomainProblem<2> bidomain_problem( &cell_factory );
bidomain_problem.Initialise();
bidomain_problem.Solve();
The results can now be visualised - the effect of the fibres changing direction at x=0.05 on the wave should be very clear.
Please see Visualisation Guides for details of using Meshalyzer/Cmgui/Paraview.
We described in the previous tutorial how to access the latest voltage vector using
ReplicatableVector
, here we illustrate how to access the voltage values using the
DistributedVector
class, which can be used to only iterate over the values
of the voltage owned by that process (for parallel simulations).
DistributedVector dist_bidomain_voltage = bidomain_problem.GetSolutionDistributedVector();
DistributedVector::Stripe bidomain_voltage(dist_bidomain_voltage, 0);
DistributedVector::Stripe extracellular_potential(dist_bidomain_voltage, 1);
A loop over all the components owned by this process..
for (DistributedVector::Iterator index = dist_bidomain_voltage.Begin();
index != dist_bidomain_voltage.End();
++index)
{
.. and a simple test, that the ’last’ node was stimulated:
if (index.Global==bidomain_problem.rGetMesh().GetNumNodes()-1) // ie if the last node
{
TS_ASSERT_LESS_THAN(0, bidomain_voltage[index]);
}
}
}
};
Full code
#include <cxxtest/TestSuite.h>
#include "BidomainProblem.hpp"
#include "LuoRudy1991.hpp"
#include "PetscSetupAndFinalize.hpp"
#include "PlaneStimulusCellFactory.hpp"
class TestAnotherBidomainTutorial : public CxxTest::TestSuite
{
public:
void TestFibreSimulation()
{
HeartConfig::Instance()->Reset();
PlaneStimulusCellFactory<CellLuoRudy1991FromCellML,2> cell_factory(-2000000);
HeartConfig::Instance()->SetSimulationDuration(5.0); //ms
HeartConfig::Instance()->SetOutputDirectory("BidomainFibresTutorial");
HeartConfig::Instance()->SetOutputFilenamePrefix("results");
HeartConfig::Instance()->SetMeshFileName("mesh/test/data/2D_0_to_1mm_800_elements", cp::media_type::Orthotropic);
// TetrahedralMesh<2,2> mesh;
// TrianglesMeshReader<2,2> reader("mesh/test/data/2D_0_to_1mm_800_elements");
// mesh.ConstructFromMeshReader(reader);
// for (unsigned i=0; i<mesh.GetNumElements(); i++)
// {
// double x = mesh.GetElement(i)->CalculateCentroid()[0];
// if (x<0.05)
// {
// std::cout << "1 0 0 1\n";
// }
// else
// {
// std::cout << "0.707106781 0.707106781 -0.707106781 0.707106781\n";
// }
// }
// return;
double scale = 2;
HeartConfig::Instance()->SetIntracellularConductivities(Create_c_vector(1.75*scale, 0.19*scale));
HeartConfig::Instance()->SetExtracellularConductivities(Create_c_vector(7.0*scale, 2.4*scale));
//HeartConfig::Instance()->SetVisualizeWithMeshalyzer(true);
HeartConfig::Instance()->SetVisualizeWithCmgui(true);
//HeartConfig::Instance()->SetVisualizeWithVtk(true);
//HeartConfig::Instance()->SetVisualizeWithParallelVtk(true);
BidomainProblem<2> bidomain_problem( &cell_factory );
bidomain_problem.Initialise();
bidomain_problem.Solve();
DistributedVector dist_bidomain_voltage = bidomain_problem.GetSolutionDistributedVector();
DistributedVector::Stripe bidomain_voltage(dist_bidomain_voltage, 0);
DistributedVector::Stripe extracellular_potential(dist_bidomain_voltage, 1);
for (DistributedVector::Iterator index = dist_bidomain_voltage.Begin();
index != dist_bidomain_voltage.End();
++index)
{
if (index.Global==bidomain_problem.rGetMesh().GetNumNodes()-1) // ie if the last node
{
TS_ASSERT_LESS_THAN(0, bidomain_voltage[index]);
}
}
}
};