This tutorial is automatically generated from the file heart/test/tutorials/TestAnotherCardiacElectroMechanicsTutorial.hpp at revision 5e8c8d7218a9/git_repo. Note that the code is given in full at the bottom of the page.
Cardiac Electro-mechanical Problems (cont.)
It is worth running this test suite with build=GccOpt_ndebug
The same includes as the previous tutorial
#include <cxxtest/TestSuite.h>
#include "PlaneStimulusCellFactory.hpp"
#include "PetscSetupAndFinalize.hpp"
#include "CardiacElectroMechProbRegularGeom.hpp"
#include "CardiacElectroMechanicsProblem.hpp"
#include "LuoRudy1991.hpp"
#include "NonlinearElasticityTools.hpp"
#include "NobleVargheseKohlNoble1998WithSac.hpp"
#include "CompressibleMooneyRivlinMaterialLaw.hpp"
#include "NobleVargheseKohlNoble1998WithSac.hpp"
#include "Hdf5ToMeshalyzerConverter.hpp"
#include "ZeroStimulusCellFactory.hpp"
#include "FileFinder.hpp"
A cell factory used in one of the tests
class PointStimulus2dCellFactory : public AbstractCardiacCellFactory<2> { private: boost::shared_ptr<SimpleStimulus> mpStimulus; public: PointStimulus2dCellFactory() : AbstractCardiacCellFactory<2>(), mpStimulus(new SimpleStimulus(-5e5, 0.5)) { } AbstractCardiacCell* CreateCardiacCellForTissueNode(Node<2>* pNode) { double x = pNode->rGetLocation()[0]; double y = pNode->rGetLocation()[1]; if (fabs(x)<0.02+1e-6 && y<-0.4+1e-6) // stimulating small region { return new CellLuoRudy1991FromCellML(mpSolver, mpStimulus); } else { return new CellLuoRudy1991FromCellML(mpSolver, mpZeroStimulus); } } }; class TestAnotherCardiacElectroMechanicsTutorial : public CxxTest::TestSuite { public:
Mechano-electric feedback, and alternative boundary conditions
Let us now run a simulation with mechano-electric feedback (MEF), and with different boundary conditions.
void TestWithMef() {
If we want to use MEF, where the stretch (in the fibre-direction) couples back to the cell model and is used in stretch-activated channels (SACs), we can't just let Chaste convert from cellml to C++ code as usual (see electro-physiology tutorials on how cell model files are autogenerated from CellML during compilation), since these files don't use stretch and don't have SACs. We have to use chaste_codegen to create a cell model class for us, rename and save it, and manually add the SAC current.
There is one example of this already in the code-base, which we will use it the following simulation. It is the Noble 98 model, with a SAC added that depends linearly on stretches (>1). It is found in the file !NobleVargheseKohlNoble1998WithSac.hpp, and is called CML_noble_varghese_kohl_noble_1998_basic_with_sac.
To add a SAC current to (or otherwise alter) your favourite cell model, you have to do the following. Auto-generate the C++ code, by running the following on the cellml file: ./python/ConvertCellModel.py heart/src/odes/cellml/LuoRudy1991.cellml (see ChasteGuides/CodeGenerationFromCellML#Usingthehelperscript if you want further documentation on this script).
Copy and rename the resultant .hpp and .cpp files (which can be found in the same folder as the input cellml). For example, rename everything to LuoRudy1991WithSac. Then alter the class to overload the method AbstractCardiacCell::SetStretch(double stretch) to store the stretch, and then implement the SAC in the GetIIonic() method. CML_noble_varghese_kohl_noble_1998_basic_with_sac provides an example of the changes that need to be made.
Let us create a cell factory returning these Noble98 SAC cells, but with no stimulus - the SAC switching on will lead be to activation.
ZeroStimulusCellFactory<CML_noble_varghese_kohl_noble_1998_basic_with_sac, 2> cell_factory;
Construct two meshes are before, in 2D
TetrahedralMesh<2,2> electrics_mesh; electrics_mesh.ConstructRegularSlabMesh(0.01/*stepsize*/, 0.1/*length*/, 0.1/*width*/, 0.1/*depth*/); QuadraticMesh<2> mechanics_mesh; mechanics_mesh.ConstructRegularSlabMesh(0.02, 0.1, 0.1, 0.1 /*as above with a different stepsize*/);
Collect the fixed nodes. This time we directly specify the new locations. We say the nodes on X=0 are to be fixed, setting the deformed x=0, but leaving y to be free (sliding boundary conditions). This functionality is described in more detail in the solid mechanics tutorials.
std::vector<unsigned> fixed_nodes; std::vector<c_vector<double,2> > fixed_node_locations; fixed_nodes.push_back(0); fixed_node_locations.push_back(zero_vector<double>(2)); for (unsigned i=1; i<mechanics_mesh.GetNumNodes(); i++) { double X = mechanics_mesh.GetNode(i)->rGetLocation()[0]; if (fabs(X) < 1e-6) // ie, if X==0 { c_vector<double,2> new_position; new_position(0) = 0.0; new_position(1) = ElectroMechanicsProblemDefinition<2>::FREE; fixed_nodes.push_back(i); fixed_node_locations.push_back(new_position); } }
Now specify tractions on the top and bottom surfaces. For full descriptions of how to apply tractions see the solid mechanics tutorials. Here, we collect the boundary elements on the bottom and top surfaces, and apply inward tractions - this will have the effect of stretching the tissue in the X-direction.
std::vector<BoundaryElement<1,2>*> boundary_elems; std::vector<c_vector<double,2> > tractions; c_vector<double,2> traction; for (TetrahedralMesh<2,2>::BoundaryElementIterator iter = mechanics_mesh.GetBoundaryElementIteratorBegin(); iter != mechanics_mesh.GetBoundaryElementIteratorEnd(); ++iter) { if (fabs((*iter)->CalculateCentroid()[1] - 0.0) < 1e-6) // if Y=0 { BoundaryElement<1,2>* p_element = *iter; boundary_elems.push_back(p_element); traction(0) = 0.0; // kPa, since the contraction model and material law use kPa for stiffnesses traction(1) = 1.0; // kPa, since the contraction model and material law use kPa for stiffnesses tractions.push_back(traction); } if (fabs((*iter)->CalculateCentroid()[1] - 0.1) < 1e-6) // if Y=0.1 { BoundaryElement<1,2>* p_element = *iter; boundary_elems.push_back(p_element); traction(0) = 0.0; traction(1) = -1.0; tractions.push_back(traction); } }
Now set up the problem. We will use a compressible approach.
ElectroMechanicsProblemDefinition<2> problem_defn(mechanics_mesh); problem_defn.SetContractionModel(KERCHOFFS2003,0.01/*contraction model ODE timestep*/); problem_defn.SetUseDefaultCardiacMaterialLaw(INCOMPRESSIBLE); problem_defn.SetMechanicsSolveTimestep(1.0);
Set the fixed node and traction info.
problem_defn.SetFixedNodes(fixed_nodes, fixed_node_locations); problem_defn.SetTractionBoundaryConditions(boundary_elems, tractions);
Now say that the deformation should affect the electro-physiology
problem_defn.SetDeformationAffectsElectrophysiology(false /*deformation affects conductivity*/, true /*deformation affects cell models*/);
Set the end time, create the problem, and solve
HeartConfig::Instance()->SetSimulationDuration(50.0); CardiacElectroMechanicsProblem<2,1> problem(INCOMPRESSIBLE, MONODOMAIN, &electrics_mesh, &mechanics_mesh, &cell_factory, &problem_defn, "TestCardiacElectroMechanicsWithMef"); problem.Solve();
Nothing exciting happens in the simulation as it is currently written. To get some interesting occurring, alter the SAC conductance in the cell model from 0.035 to 0.35 (micro-Siemens). (look for the line const double g_sac = 0.035 in NobleVargheseKohlNoble1998WithSac.hpp).
Rerun and visualise as usual, using Cmgui. By visualising the voltage on the deforming mesh, you can see that the voltage gradually increases due to the SAC, since the tissue is stretched, until the threshold is reached and activation occurs.
For MEF simulations, we may want to visualise the electrical results on the electrics mesh using Meshalyzer, for example to more easily visualise action potentials. This isn't (and currently can't be) created by CardiacElectroMechanicsProblem. We can use a converter as follows to post-process:
FileFinder test_output_folder("TestCardiacElectroMechanicsWithMef/electrics", RelativeTo::ChasteTestOutput); Hdf5ToMeshalyzerConverter<2,2> converter(test_output_folder, "voltage", &electrics_mesh, false, HeartConfig::Instance()->GetVisualizerOutputPrecision());
Some other notes. If you want to apply time-dependent traction boundary conditions, this is possible by specifying the traction in functional form - see solid mechanics tutorials. Similarly, more natural 'pressure acting on the deformed body' boundary conditions are possible - see below tutorial.
Robustness: Sometimes the nonlinear solver doesn't converge, and will give an error. This can be due to either a non-physical (or not very physical) situation, or just because the current guess is quite far from the solution and the solver can't find the solution (this can occur in nonlinear elasticity problems if the loading is large, for example). Current work in progress is on making the solver more robust, and also on parallelising the solver. One option when a solve fails is to decrease the mechanics timestep.
Ignore the following, it is just to check nothing has changed.
Hdf5DataReader reader("TestCardiacElectroMechanicsWithMef/electrics", "voltage"); unsigned num_timesteps = reader.GetUnlimitedDimensionValues().size(); Vec voltage = PetscTools::CreateVec(electrics_mesh.GetNumNodes()); reader.GetVariableOverNodes(voltage, "V", num_timesteps-1); ReplicatableVector voltage_repl(voltage); for (unsigned i=0; i<voltage_repl.GetSize(); i++) { TS_ASSERT_DELTA(voltage_repl[i], -81.9080, 1e-3); } PetscTools::Destroy(voltage); }
Internal pressures
Next, we run a simulation on a 2d annulus, with an internal pressure applied.
void TestAnnulusWithInternalPressure() {
The following should require little explanation now
TetrahedralMesh<2,2> electrics_mesh; QuadraticMesh<2> mechanics_mesh; // could (should?) use finer electrics mesh, but keeping electrics simulation time down TrianglesMeshReader<2,2> reader1("mesh/test/data/annuli/circular_annulus_960_elements"); electrics_mesh.ConstructFromMeshReader(reader1); TrianglesMeshReader<2,2> reader2("mesh/test/data/annuli/circular_annulus_960_elements_quad",2 /*quadratic elements*/); mechanics_mesh.ConstructFromMeshReader(reader2); PointStimulus2dCellFactory cell_factory; std::vector<unsigned> fixed_nodes; std::vector<c_vector<double,2> > fixed_node_locations; for (unsigned i=0; i<mechanics_mesh.GetNumNodes(); i++) { double x = mechanics_mesh.GetNode(i)->rGetLocation()[0]; double y = mechanics_mesh.GetNode(i)->rGetLocation()[1]; if (fabs(x)<1e-6 && fabs(y+0.5)<1e-6) // fixed point (0.0,-0.5) at bottom of mesh { fixed_nodes.push_back(i); c_vector<double,2> new_position; new_position(0) = x; new_position(1) = y; fixed_node_locations.push_back(new_position); } if (fabs(x)<1e-6 && fabs(y-0.5)<1e-6) // constrained point (0.0,0.5) at top of mesh { fixed_nodes.push_back(i); c_vector<double,2> new_position; new_position(0) = x; new_position(1) = ElectroMechanicsProblemDefinition<2>::FREE; fixed_node_locations.push_back(new_position); } }
Increase this end time to see more contraction
HeartConfig::Instance()->SetSimulationDuration(30.0); ElectroMechanicsProblemDefinition<2> problem_defn(mechanics_mesh); problem_defn.SetContractionModel(KERCHOFFS2003,0.1); problem_defn.SetUseDefaultCardiacMaterialLaw(COMPRESSIBLE); //problem_defn.SetZeroDisplacementNodes(fixed_nodes); problem_defn.SetFixedNodes(fixed_nodes, fixed_node_locations); problem_defn.SetMechanicsSolveTimestep(1.0); FileFinder finder("heart/test/data/fibre_tests/circular_annulus_960_elements.ortho",RelativeTo::ChasteSourceRoot); problem_defn.SetVariableFibreSheetDirectionsFile(finder, false);
The elasticity solvers have two nonlinear solvers implemented, one hand-coded and one which uses PETSc's SNES solver. The latter is not the default but can be more robust (and will probably be the default in later versions). This is how it can be used. (This option can also be called if the compiled binary is run from the command line (see ChasteGuides/RunningBinariesFromCommandLine) using the option "-mech_use_snes").
problem_defn.SetSolveUsingSnes();
Now let us collect all the boundary elements on the inner (endocardial) surface. The following uses knowledge about the geometry - the inner surface is r=0.3, the outer is r=0.5.
std::vector<BoundaryElement<1,2>*> boundary_elems; for (TetrahedralMesh<2,2>::BoundaryElementIterator iter = mechanics_mesh.GetBoundaryElementIteratorBegin(); iter != mechanics_mesh.GetBoundaryElementIteratorEnd(); ++iter) { ChastePoint<2> centroid = (*iter)->CalculateCentroid(); double r = sqrt( centroid[0]*centroid[0] + centroid[1]*centroid[1] ); if (r < 0.4) { BoundaryElement<1,2>* p_element = *iter; boundary_elems.push_back(p_element); } }
This is how to set the pressure to be applied to these boundary elements. The negative sign implies inward pressure.
problem_defn.SetApplyNormalPressureOnDeformedSurface(boundary_elems, -1.0 /*1 KPa is about 8mmHg*/);
The solver computes the equilibrium solution (given the pressure loading) before the first timestep. As there is a big deformation from the undeformed state to this loaded state, the nonlinear solver may not converge. The following increments the loading (solves with p=-1/3, then p=-2/3, then p=-1), which allows convergence to occur.
problem_defn.SetNumIncrementsForInitialDeformation(3); CardiacElectroMechanicsProblem<2,1> problem(COMPRESSIBLE, MONODOMAIN, &electrics_mesh, &mechanics_mesh, &cell_factory, &problem_defn, "TestAnnulusWithInternalPressure");
If we want stresses and strains output, we can do the following. The deformation gradients and 2nd PK stresses for each element will be written at the requested times.
problem.SetOutputDeformationGradientsAndStress(10.0 /*how often (in ms) to write - should be a multiple of mechanics timestep*/);
Since this test involves a large deformation at t=0, several Newton iterations are required. To see how the nonlinear solve is progressing, you can run from the binary from the command line with the command line argument "-mech_verbose".
problem.Solve();
Visualise using cmgui, and note the different shapes at t=-1 (undeformed) and t=0 (loaded)
Note: if you want to have a time-dependent pressure, you can replace the second parameter (the pressure) in SetApplyNormalPressureOnDeformedSurface() with a function pointer (the name of a function) which returns the pressure as a function of time.
}
More examples: For a 3d ellipsoid geometry test, see heart/test/mechanics/TestCardiacElectroMechanicsOnEllipsoid.hpp
};
Code
The full code is given below
File name TestAnotherCardiacElectroMechanicsTutorial.hpp
#include <cxxtest/TestSuite.h> #include "PlaneStimulusCellFactory.hpp" #include "PetscSetupAndFinalize.hpp" #include "CardiacElectroMechProbRegularGeom.hpp" #include "CardiacElectroMechanicsProblem.hpp" #include "LuoRudy1991.hpp" #include "NonlinearElasticityTools.hpp" #include "NobleVargheseKohlNoble1998WithSac.hpp" #include "CompressibleMooneyRivlinMaterialLaw.hpp" #include "NobleVargheseKohlNoble1998WithSac.hpp" #include "Hdf5ToMeshalyzerConverter.hpp" #include "ZeroStimulusCellFactory.hpp" #include "FileFinder.hpp" class PointStimulus2dCellFactory : public AbstractCardiacCellFactory<2> { private: boost::shared_ptr<SimpleStimulus> mpStimulus; public: PointStimulus2dCellFactory() : AbstractCardiacCellFactory<2>(), mpStimulus(new SimpleStimulus(-5e5, 0.5)) { } AbstractCardiacCell* CreateCardiacCellForTissueNode(Node<2>* pNode) { double x = pNode->rGetLocation()[0]; double y = pNode->rGetLocation()[1]; if (fabs(x)<0.02+1e-6 && y<-0.4+1e-6) // stimulating small region { return new CellLuoRudy1991FromCellML(mpSolver, mpStimulus); } else { return new CellLuoRudy1991FromCellML(mpSolver, mpZeroStimulus); } } }; class TestAnotherCardiacElectroMechanicsTutorial : public CxxTest::TestSuite { public: void TestWithMef() { ZeroStimulusCellFactory<CML_noble_varghese_kohl_noble_1998_basic_with_sac, 2> cell_factory; TetrahedralMesh<2,2> electrics_mesh; electrics_mesh.ConstructRegularSlabMesh(0.01/*stepsize*/, 0.1/*length*/, 0.1/*width*/, 0.1/*depth*/); QuadraticMesh<2> mechanics_mesh; mechanics_mesh.ConstructRegularSlabMesh(0.02, 0.1, 0.1, 0.1 /*as above with a different stepsize*/); std::vector<unsigned> fixed_nodes; std::vector<c_vector<double,2> > fixed_node_locations; fixed_nodes.push_back(0); fixed_node_locations.push_back(zero_vector<double>(2)); for (unsigned i=1; i<mechanics_mesh.GetNumNodes(); i++) { double X = mechanics_mesh.GetNode(i)->rGetLocation()[0]; if (fabs(X) < 1e-6) // ie, if X==0 { c_vector<double,2> new_position; new_position(0) = 0.0; new_position(1) = ElectroMechanicsProblemDefinition<2>::FREE; fixed_nodes.push_back(i); fixed_node_locations.push_back(new_position); } } std::vector<BoundaryElement<1,2>*> boundary_elems; std::vector<c_vector<double,2> > tractions; c_vector<double,2> traction; for (TetrahedralMesh<2,2>::BoundaryElementIterator iter = mechanics_mesh.GetBoundaryElementIteratorBegin(); iter != mechanics_mesh.GetBoundaryElementIteratorEnd(); ++iter) { if (fabs((*iter)->CalculateCentroid()[1] - 0.0) < 1e-6) // if Y=0 { BoundaryElement<1,2>* p_element = *iter; boundary_elems.push_back(p_element); traction(0) = 0.0; // kPa, since the contraction model and material law use kPa for stiffnesses traction(1) = 1.0; // kPa, since the contraction model and material law use kPa for stiffnesses tractions.push_back(traction); } if (fabs((*iter)->CalculateCentroid()[1] - 0.1) < 1e-6) // if Y=0.1 { BoundaryElement<1,2>* p_element = *iter; boundary_elems.push_back(p_element); traction(0) = 0.0; traction(1) = -1.0; tractions.push_back(traction); } } ElectroMechanicsProblemDefinition<2> problem_defn(mechanics_mesh); problem_defn.SetContractionModel(KERCHOFFS2003,0.01/*contraction model ODE timestep*/); problem_defn.SetUseDefaultCardiacMaterialLaw(INCOMPRESSIBLE); problem_defn.SetMechanicsSolveTimestep(1.0); problem_defn.SetFixedNodes(fixed_nodes, fixed_node_locations); problem_defn.SetTractionBoundaryConditions(boundary_elems, tractions); problem_defn.SetDeformationAffectsElectrophysiology(false /*deformation affects conductivity*/, true /*deformation affects cell models*/); HeartConfig::Instance()->SetSimulationDuration(50.0); CardiacElectroMechanicsProblem<2,1> problem(INCOMPRESSIBLE, MONODOMAIN, &electrics_mesh, &mechanics_mesh, &cell_factory, &problem_defn, "TestCardiacElectroMechanicsWithMef"); problem.Solve(); FileFinder test_output_folder("TestCardiacElectroMechanicsWithMef/electrics", RelativeTo::ChasteTestOutput); Hdf5ToMeshalyzerConverter<2,2> converter(test_output_folder, "voltage", &electrics_mesh, false, HeartConfig::Instance()->GetVisualizerOutputPrecision()); Hdf5DataReader reader("TestCardiacElectroMechanicsWithMef/electrics", "voltage"); unsigned num_timesteps = reader.GetUnlimitedDimensionValues().size(); Vec voltage = PetscTools::CreateVec(electrics_mesh.GetNumNodes()); reader.GetVariableOverNodes(voltage, "V", num_timesteps-1); ReplicatableVector voltage_repl(voltage); for (unsigned i=0; i<voltage_repl.GetSize(); i++) { TS_ASSERT_DELTA(voltage_repl[i], -81.9080, 1e-3); } PetscTools::Destroy(voltage); } void TestAnnulusWithInternalPressure() { TetrahedralMesh<2,2> electrics_mesh; QuadraticMesh<2> mechanics_mesh; // could (should?) use finer electrics mesh, but keeping electrics simulation time down TrianglesMeshReader<2,2> reader1("mesh/test/data/annuli/circular_annulus_960_elements"); electrics_mesh.ConstructFromMeshReader(reader1); TrianglesMeshReader<2,2> reader2("mesh/test/data/annuli/circular_annulus_960_elements_quad",2 /*quadratic elements*/); mechanics_mesh.ConstructFromMeshReader(reader2); PointStimulus2dCellFactory cell_factory; std::vector<unsigned> fixed_nodes; std::vector<c_vector<double,2> > fixed_node_locations; for (unsigned i=0; i<mechanics_mesh.GetNumNodes(); i++) { double x = mechanics_mesh.GetNode(i)->rGetLocation()[0]; double y = mechanics_mesh.GetNode(i)->rGetLocation()[1]; if (fabs(x)<1e-6 && fabs(y+0.5)<1e-6) // fixed point (0.0,-0.5) at bottom of mesh { fixed_nodes.push_back(i); c_vector<double,2> new_position; new_position(0) = x; new_position(1) = y; fixed_node_locations.push_back(new_position); } if (fabs(x)<1e-6 && fabs(y-0.5)<1e-6) // constrained point (0.0,0.5) at top of mesh { fixed_nodes.push_back(i); c_vector<double,2> new_position; new_position(0) = x; new_position(1) = ElectroMechanicsProblemDefinition<2>::FREE; fixed_node_locations.push_back(new_position); } } HeartConfig::Instance()->SetSimulationDuration(30.0); ElectroMechanicsProblemDefinition<2> problem_defn(mechanics_mesh); problem_defn.SetContractionModel(KERCHOFFS2003,0.1); problem_defn.SetUseDefaultCardiacMaterialLaw(COMPRESSIBLE); //problem_defn.SetZeroDisplacementNodes(fixed_nodes); problem_defn.SetFixedNodes(fixed_nodes, fixed_node_locations); problem_defn.SetMechanicsSolveTimestep(1.0); FileFinder finder("heart/test/data/fibre_tests/circular_annulus_960_elements.ortho",RelativeTo::ChasteSourceRoot); problem_defn.SetVariableFibreSheetDirectionsFile(finder, false); problem_defn.SetSolveUsingSnes(); std::vector<BoundaryElement<1,2>*> boundary_elems; for (TetrahedralMesh<2,2>::BoundaryElementIterator iter = mechanics_mesh.GetBoundaryElementIteratorBegin(); iter != mechanics_mesh.GetBoundaryElementIteratorEnd(); ++iter) { ChastePoint<2> centroid = (*iter)->CalculateCentroid(); double r = sqrt( centroid[0]*centroid[0] + centroid[1]*centroid[1] ); if (r < 0.4) { BoundaryElement<1,2>* p_element = *iter; boundary_elems.push_back(p_element); } } problem_defn.SetApplyNormalPressureOnDeformedSurface(boundary_elems, -1.0 /*1 KPa is about 8mmHg*/); problem_defn.SetNumIncrementsForInitialDeformation(3); CardiacElectroMechanicsProblem<2,1> problem(COMPRESSIBLE, MONODOMAIN, &electrics_mesh, &mechanics_mesh, &cell_factory, &problem_defn, "TestAnnulusWithInternalPressure"); problem.SetOutputDeformationGradientsAndStress(10.0 /*how often (in ms) to write - should be a multiple of mechanics timestep*/); problem.Solve(); } };