37#ifndef ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
38#define ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
44#include <boost/serialization/base_object.hpp>
47#include "AbstractCardiacCell.hpp"
50#include "TimeStepper.hpp"
66template<
unsigned SIZE>
78 template<
class Archive>
79 void serialize(Archive & archive,
const unsigned int version)
82 archive & boost::serialization::base_object<AbstractCardiacCell>(*
this);
102 unsigned numberOfStateVariables,
103 unsigned voltageIndex,
104 boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus);
116 virtual void ComputeResidual(
double time,
const double rCurrentGuess[SIZE],
double rResidual[SIZE])=0;
125 virtual void ComputeJacobian(
double time,
const double rCurrentGuess[SIZE],
double rJacobian[SIZE][SIZE])=0;
207template <
unsigned SIZE>
209 unsigned numberOfStateVariables,
210 unsigned voltageIndex,
211 boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
213 numberOfStateVariables,
215 pIntracellularStimulus)
218template <
unsigned SIZE>
222template <
unsigned SIZE>
235 double _n_steps = (tEnd - tStart) / tSamp;
236 const unsigned n_steps = (
unsigned) floor(_n_steps+0.5);
237 assert(fabs(tStart+n_steps*tSamp - tEnd) < 1e-12);
238 const unsigned n_small_steps = (
unsigned) floor(tSamp/mDt+0.5);
239 assert(fabs(mDt*n_small_steps - tSamp) < 1e-12);
249 double curr_time = tStart;
250 for (
unsigned i=0; i<n_steps; i++)
252 for (
unsigned j=0; j<n_small_steps; j++)
254 curr_time = tStart + i*tSamp + j*mDt;
257 UpdateTransmembranePotential(curr_time);
260 ComputeOneStepExceptVoltage(curr_time);
263 VerifyStateVariables();
268 solutions.
rGetTimes().push_back(curr_time+mDt);
274template <
unsigned SIZE>
281 unsigned n_steps = (
unsigned)((tEnd - tStart) / mDt + 0.5);
282 assert(fabs(tStart + n_steps*mDt - tEnd) < 1e-12);
286 for (
unsigned i=0; i<n_steps; i++)
288 curr_time = tStart + i*mDt;
291 ComputeOneStepExceptVoltage(curr_time);
295 VerifyStateVariables();
300template<
unsigned SIZE>
307 double time = stepper.
GetTime();
310 UpdateTransmembranePotential(time);
313 ComputeOneStepExceptVoltage(time);
316 VerifyStateVariables();
334 friend class boost::serialization::access;
341 template<
class Archive>
342 void serialize(Archive & archive,
const unsigned int version)
345 archive & boost::serialization::base_object<AbstractCardiacCell>(*
this);
358 unsigned voltageIndex,
359 boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
361 numberOfStateVariables,
363 pIntracellularStimulus)
394 double _n_steps = (tEnd - tStart) / tSamp;
395 const unsigned n_steps = (
unsigned) floor(_n_steps+0.5);
396 assert(fabs(tStart+n_steps*tSamp - tEnd) < 1e-12);
397 const unsigned n_small_steps = (
unsigned) floor(tSamp/
mDt+0.5);
398 assert(fabs(
mDt*n_small_steps - tSamp) < 1e-12);
408 double curr_time = tStart;
409 for (
unsigned i=0; i<n_steps; i++)
411 for (
unsigned j=0; j<n_small_steps; j++)
413 curr_time = tStart + i*tSamp + j*
mDt;
448 unsigned n_steps = (
unsigned)((tEnd - tStart) /
mDt + 0.5);
449 assert(fabs(tStart + n_steps*
mDt - tEnd) < 1e-12);
453 for (
unsigned i=0; i<n_steps; i++)
455 curr_time = tStart + i*
mDt;
480 double time = stepper.
GetTime();
534#include "OutputFileHandler.hpp"
535#include <boost/foreach.hpp>
536template<
unsigned SIZE>
537void DumpJacobianToFile(
double time,
const double rCurrentGuess[SIZE],
double rJacobian[SIZE][SIZE],
538 const std::vector<double>& rY)
541 out_stream p_file = handler.
OpenOutputFile(
"J.txt", std::ios::app);
542 (*p_file) <<
"At " << time <<
" " << SIZE << std::endl;
544 BOOST_FOREACH(
double y, rY)
546 (*p_file) <<
" " << y;
548 (*p_file) << std::endl;
549 (*p_file) <<
"rCurrentGuess";
550 for (
unsigned i=0; i<SIZE; i++)
552 (*p_file) <<
" " << rCurrentGuess[i];
554 (*p_file) << std::endl;
555 (*p_file) <<
"rJacobian";
556 for (
unsigned i=0; i<SIZE; i++)
558 for (
unsigned j=0; j<SIZE; j++)
560 (*p_file) <<
" " << rJacobian[i][j];
563 (*p_file) << std::endl;
#define TEMPLATED_CLASS_IS_ABSTRACT_1_UNSIGNED(T)
AbstractBackwardEulerCardiacCell(unsigned numberOfStateVariables, unsigned voltageIndex, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
void ComputeExceptVoltage(double tStart, double tEnd)
OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)
virtual ~AbstractBackwardEulerCardiacCell()
void SolveAndUpdateState(double tStart, double tEnd)
void serialize(Archive &archive, const unsigned int version)
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
virtual void ComputeOneStepExceptVoltage(double tStart)=0
virtual void UpdateTransmembranePotential(double time)=0
void ComputeExceptVoltage(double tStart, double tEnd)
virtual void ComputeJacobian(double time, const double rCurrentGuess[SIZE], double rJacobian[SIZE][SIZE])=0
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
virtual ~AbstractBackwardEulerCardiacCell()
virtual void ComputeOneStepExceptVoltage(double tStart)=0
void serialize(Archive &archive, const unsigned int version)
virtual void ComputeResidual(double time, const double rCurrentGuess[SIZE], double rResidual[SIZE])=0
void SolveAndUpdateState(double tStart, double tEnd)
friend class boost::serialization::access
OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)
virtual void UpdateTransmembranePotential(double time)=0
AbstractBackwardEulerCardiacCell(unsigned numberOfStateVariables, unsigned voltageIndex, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
virtual void VerifyStateVariables()
std::vector< double > & rGetStateVariables()
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
void SetNumberOfTimeSteps(unsigned numTimeSteps)
std::vector< std::vector< double > > & rGetSolutions()
std::vector< double > & rGetTimes()
void SetOdeSystemInformation(boost::shared_ptr< const AbstractOdeSystemInformation > pOdeSystemInfo)
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
void AdvanceOneTimeStep()