00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #ifndef ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
00031 #define ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
00032
00033 #include "AbstractCardiacCell.hpp"
00034
00035
00036 #include <cassert>
00037 #include <cmath>
00038
00053 template<unsigned SIZE>
00054 class AbstractBackwardEulerCardiacCell : public AbstractCardiacCell
00055 {
00056 public:
00057
00075 AbstractBackwardEulerCardiacCell(
00076 unsigned numberOfStateVariables,
00077 unsigned voltageIndex,
00078 AbstractStimulusFunction* intracellularStimulus);
00079
00080 virtual ~AbstractBackwardEulerCardiacCell();
00081
00088 virtual void ComputeResidual(const double rCurrentGuess[SIZE], double rResidual[SIZE])=0;
00089
00096 virtual void ComputeJacobian(const double rCurrentGuess[SIZE], double rJacobian[SIZE][SIZE])=0;
00097
00107 virtual OdeSolution Compute(double tStart, double tEnd);
00108
00115 virtual void ComputeExceptVoltage(double tStart, double tEnd);
00116
00117 private:
00118 #define COVERAGE_IGNORE
00119
00122 void EvaluateYDerivatives(double time, const std::vector<double> &rY, std::vector<double> &rDY)
00123 {
00124 NEVER_REACHED;
00125 }
00126 #undef COVERAGE_IGNORE
00127
00128 protected:
00135 virtual void ComputeOneStepExceptVoltage(double tStart)=0;
00136
00142 virtual void UpdateTransmembranePotential(double time)=0;
00143 };
00144
00145 template <unsigned SIZE>
00146 AbstractBackwardEulerCardiacCell<SIZE>::AbstractBackwardEulerCardiacCell(
00147 unsigned numberOfStateVariables,
00148 unsigned voltageIndex,
00149 AbstractStimulusFunction* intracellularStimulus)
00150 : AbstractCardiacCell(NULL,
00151 numberOfStateVariables,
00152 voltageIndex,
00153 intracellularStimulus)
00154 {}
00155
00156 template <unsigned SIZE>
00157 AbstractBackwardEulerCardiacCell<SIZE>::~AbstractBackwardEulerCardiacCell()
00158 {}
00159
00160 template <unsigned SIZE>
00161 OdeSolution AbstractBackwardEulerCardiacCell<SIZE>::Compute(double tStart, double tEnd)
00162 {
00163
00164
00165
00166
00167
00168 double _n_steps = (tEnd - tStart) / mDt;
00169 unsigned n_steps = (unsigned) round(_n_steps);
00170 assert(fabs(tStart+n_steps*mDt - tEnd) < 1e-12);
00171
00172
00173 OdeSolution solutions;
00174 solutions.SetNumberOfTimeSteps(n_steps);
00175 solutions.rGetSolutions().push_back(rGetStateVariables());
00176 solutions.rGetTimes().push_back(tStart);
00177
00178
00179 double curr_time;
00180 for (unsigned i=0; i<n_steps; i++)
00181 {
00182 curr_time = tStart + i*mDt;
00183
00184
00185 UpdateTransmembranePotential(curr_time);
00186
00187
00188 ComputeOneStepExceptVoltage(curr_time);
00189
00190
00191 solutions.rGetSolutions().push_back(rGetStateVariables());
00192 solutions.rGetTimes().push_back(curr_time+mDt);
00193
00194
00195 VerifyStateVariables();
00196 }
00197
00198 return solutions;
00199 }
00200
00201 template <unsigned SIZE>
00202 void AbstractBackwardEulerCardiacCell<SIZE>::ComputeExceptVoltage(double tStart, double tEnd)
00203 {
00204
00205
00206
00207 double _n_steps = (tEnd - tStart) / mDt;
00208 unsigned n_steps = (unsigned) round(_n_steps);
00209 assert(fabs(tStart+n_steps*mDt - tEnd) < 1e-12);
00210
00211
00212 double curr_time;
00213 for (unsigned i=0; i<n_steps; i++)
00214 {
00215 curr_time = tStart + i*mDt;
00216
00217
00218 ComputeOneStepExceptVoltage(curr_time);
00219
00220
00221 VerifyStateVariables();
00222 }
00223 }
00224
00225
00226 #endif