40#include "AbstractCellBasedSimulation.hpp"
41#include "CellBasedEventHandler.hpp"
43#include "ExecutableSupport.hpp"
44#include "AbstractPdeModifier.hpp"
45#include "CellDivisionLocationsWriter.hpp"
46#include "CellRemovalLocationsWriter.hpp"
48template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
50 bool deleteCellPopulationInDestructor,
54 mrCellPopulation(rCellPopulation),
55 mDeleteCellPopulationInDestructor(deleteCellPopulationInDestructor),
56 mInitialiseCells(initialiseCells),
58 mUpdateCellPopulation(true),
60 mSimulationOutputDirectory(mOutputDirectory),
63 mOutputDivisionLocations(false),
64 mOutputCellVelocities(false),
65 mSamplingTimestepMultiple(1),
66 mUpdatingTimestepMultiple(1)
83template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
86 if (mDeleteCellPopulationInDestructor)
88 delete &mrCellPopulation;
92template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
100 unsigned num_births_this_step = 0;
104 cell_iter != mrCellPopulation.End();
108 double cell_age = cell_iter->GetAge();
111 if (cell_iter->ReadyToDivide())
114 if (mrCellPopulation.IsRoomToDivide(*cell_iter))
117 unsigned parent_cell_id = cell_iter->GetCellId();
120 CellPtr p_new_cell = cell_iter->Divide();
128 if (mrCellPopulation.template HasWriter<CellDivisionLocationsWriter>())
130 c_vector<double, SPACE_DIM> cell_location = mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
132 std::stringstream division_info;
134 for (
unsigned i = 0; i < SPACE_DIM; i++)
136 division_info << cell_location[i] <<
"\t";
138 division_info <<
"\t" << cell_age <<
"\t" << parent_cell_id <<
"\t" << cell_iter->GetCellId() <<
"\t" << p_new_cell->GetCellId() <<
"\t";
140 mrCellPopulation.AddDivisionInformation(division_info.str());
144 mrCellPopulation.AddCell(p_new_cell, *cell_iter);
147 num_births_this_step++;
152 return num_births_this_step;
155template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
158 unsigned num_deaths_this_step = 0;
165 killer_iter != mCellKillers.end();
168 (*killer_iter)->CheckAndLabelCellsForApoptosisOrDeath();
171 num_deaths_this_step += mrCellPopulation.RemoveDeadCells();
173 return num_deaths_this_step;
176template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
183template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
189template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
195template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
201template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
208template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
211 mOutputDirectory = outputDirectory;
212 mSimulationOutputDirectory = mOutputDirectory;
215template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
218 return mOutputDirectory;
221template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
224 assert(samplingTimestepMultiple > 0);
225 mSamplingTimestepMultiple = samplingTimestepMultiple;
228template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
231 assert(updatingTimestepMultiple > 0);
232 mUpdatingTimestepMultiple = updatingTimestepMultiple;
235template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
238 return mrCellPopulation;
241template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
244 return mrCellPopulation;
247template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
250 mUpdateCellPopulation = updateCellPopulation;
253template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
256 return mUpdateCellPopulation;
259template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
265template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
268 mCellKillers.push_back(pCellKiller);
271template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
274 mCellKillers.clear();
277template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
280 mSimulationModifiers.push_back(pSimulationModifier);
283template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
286 return &mSimulationModifiers;
289template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
292 mTopologyUpdateSimulationModifiers.push_back(pSimulationModifier);
295template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
298 return &mTopologyUpdateSimulationModifiers;
301template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
304 std::vector<double> location;
305 for (
unsigned i=0; i<SPACE_DIM; i++)
307 location.push_back(mrCellPopulation.GetNode(rNodeIndex)->rGetLocation()[i]);
312template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
320 double current_time = p_simulation_time->
GetTime();
336 unsigned num_time_steps = (
unsigned) ((mEndTime-current_time)/mDt+0.5);
337 if (current_time > 0)
345 EXCEPTION(
"End time and number of timesteps already setup. You should not use SimulationTime::SetEndTimeAndNumberOfTimeSteps in cell-based tests.");
353 if (mOutputDirectory ==
"")
358 double time_now = p_simulation_time->
GetTime();
359 std::ostringstream time_string;
360 time_string << time_now;
362 std::string results_directory = mOutputDirectory +
"/results_from_time_" + time_string.str();
363 mSimulationOutputDirectory = results_directory;
370 if (mOutputDivisionLocations)
372 mrCellPopulation.template AddCellPopulationEventWriter<CellDivisionLocationsWriter>();
373 mrCellPopulation.template AddCellPopulationEventWriter<CellRemovalLocationsWriter>();
375 if (mOutputCellVelocities)
378 mpCellVelocitiesFile = output_file_handler2.OpenOutputFile(
"cellvelocities.dat");
381 mrCellPopulation.OpenWritersFiles(output_file_handler);
385 mpVizSetupFile = output_file_handler.OpenOutputFile(
"results.vizsetup");
388 iter != mSimulationModifiers.end();
393 *this->mpVizSetupFile <<
"PDE \n";
398 this->mrCellPopulation.SimulationSetupHook(
this);
404 iter != mSimulationModifiers.end();
407 (*iter)->SetupSolve(this->mrCellPopulation,this->mSimulationOutputDirectory);
412 iter != mTopologyUpdateSimulationModifiers.end();
415 (*iter)->SetupSolve(this->mrCellPopulation,this->mSimulationOutputDirectory);
422 LOG(1,
"Setting up cells...");
424 cell_iter != mrCellPopulation.End();
431 cell_iter->ReadyToDivide();
436 WriteVisualizerSetupFile();
440 *mpVizSetupFile << std::flush;
443 mrCellPopulation.WriteResultsToFiles(results_directory+
"/");
445 OutputSimulationSetup();
449 while (!( p_simulation_time->
IsFinished() || StoppingEventHasOccurred() ) )
451 LOG(1,
"--TIME = " << p_simulation_time->
GetTime() <<
"\n");
457 UpdateCellPopulation();
462 bool at_sampling_timestep = (p_time->
GetTimeStepsElapsed()%this->mSamplingTimestepMultiple == 0);
471 std::map<CellPtr, c_vector<double, SPACE_DIM> > old_cell_locations;
472 if (mOutputCellVelocities && at_sampling_timestep)
475 cell_iter != mrCellPopulation.End();
478 old_cell_locations[*cell_iter] = mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
485 iter != mTopologyUpdateSimulationModifiers.end();
488 (*iter)->UpdateAtEndOfTimeStep(this->mrCellPopulation);
493 UpdateCellLocationsAndTopology();
496 if (mOutputCellVelocities && at_sampling_timestep)
499 *mpCellVelocitiesFile << p_time->
GetTime() + mDt<<
"\t";
502 cell_iter != mrCellPopulation.End();
505 unsigned index = mrCellPopulation.GetLocationIndexUsingCell(*cell_iter);
506 const c_vector<double,SPACE_DIM>& position = mrCellPopulation.GetLocationOfCellCentre(*cell_iter);
508 c_vector<double, SPACE_DIM> velocity;
509 velocity = (position - old_cell_locations[*cell_iter])/mDt;
511 *mpCellVelocitiesFile << index <<
" ";
512 for (
unsigned i=0; i<SPACE_DIM; i++)
514 *mpCellVelocitiesFile << position[i] <<
" ";
517 for (
unsigned i=0; i<SPACE_DIM; i++)
519 *mpCellVelocitiesFile << velocity[i] <<
" ";
522 *mpCellVelocitiesFile <<
"\n";
526 mrCellPopulation.UpdateCellProcessLocation();
534 iter != mSimulationModifiers.end();
537 (*iter)->UpdateAtEndOfTimeStep(this->mrCellPopulation);
545 mrCellPopulation.WriteResultsToFiles(results_directory+
"/");
549 iter != mSimulationModifiers.end();
552 (*iter)->UpdateAtEndOfOutputTimeStep(this->mrCellPopulation);
558 LOG(1,
"--END TIME = " << p_simulation_time->
GetTime() <<
"\n");
565 UpdateCellPopulation();
570 iter != mSimulationModifiers.end();
573 (*iter)->UpdateAtEndOfSolve(this->mrCellPopulation);
579 mrCellPopulation.CloseWritersFiles();
581 if (mOutputCellVelocities)
583 mpCellVelocitiesFile->close();
588 *mpVizSetupFile <<
"Complete\n";
589 mpVizSetupFile->close();
596template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
602template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
607 unsigned deaths_this_step = DoCellRemoval();
608 mNumDeaths += deaths_this_step;
609 LOG(1,
"\tNum deaths = " << mNumDeaths <<
"\n");
614 unsigned births_this_step = DoCellBirth();
615 mNumBirths += births_this_step;
616 LOG(1,
"\tNum births = " << mNumBirths <<
"\n");
620 bool births_or_death_occurred = ((births_this_step>0) || (deaths_this_step>0));
625 if (mUpdateCellPopulation && (p_time->
GetTimeStepsElapsed() % mUpdatingTimestepMultiple == 0) )
627 LOG(1,
"\tUpdating cell population...");
628 mrCellPopulation.Update(births_or_death_occurred);
631 else if (births_or_death_occurred)
633 if (!mUpdateCellPopulation)
635 EXCEPTION(
"CellPopulation has had births or deaths but mUpdateCellPopulation is set to false, please set it to true.");
639 EXCEPTION(
"CellPopulation has had births or deaths but you were on a non update step, make sure your cell cycle model and killer only operate on update steps.");
649template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
652 return mOutputDivisionLocations;
655template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
658 mOutputDivisionLocations = outputDivisionLocations;
661template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
664 return mOutputCellVelocities;
667template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
670 mOutputCellVelocities = outputCellVelocities;
673template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
676 OutputFileHandler output_file_handler(this->mSimulationOutputDirectory +
"/",
false);
685 out_stream build_info_file = output_file_handler.
OpenOutputFile(
"build.info");
686 std::string build_info;
688 *build_info_file << build_info;
689 build_info_file->close();
692 out_stream parameter_file = output_file_handler.
OpenOutputFile(
"results.parameters");
695 std::string simulation_type = GetIdentifier();
697 *parameter_file <<
"<Chaste>\n";
698 *parameter_file <<
"\n\t<" << simulation_type <<
">\n";
699 OutputSimulationParameters(parameter_file);
700 *parameter_file <<
"\t</" << simulation_type <<
">\n";
701 *parameter_file <<
"\n";
704 mrCellPopulation.OutputCellPopulationInfo(parameter_file);
707 *parameter_file <<
"\n\t<CellKillers>\n";
709 iter != mCellKillers.end();
713 (*iter)->OutputCellKillerInfo(parameter_file);
715 *parameter_file <<
"\t</CellKillers>\n";
718 *parameter_file <<
"\n\t<SimulationModifiers>\n";
720 iter != mSimulationModifiers.end();
724 (*iter)->OutputSimulationModifierInfo(parameter_file);
726 *parameter_file <<
"\t</SimulationModifiers>\n";
729 OutputAdditionalSimulationSetup(parameter_file);
731 *parameter_file <<
"\n</Chaste>\n";
732 parameter_file->close();
736template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
739 *rParamsFile <<
"\t\t<Dt>" << mDt <<
"</Dt>\n";
740 *rParamsFile <<
"\t\t<EndTime>" << mEndTime <<
"</EndTime>\n";
741 *rParamsFile <<
"\t\t<UpdateCellPopulation>" << mUpdateCellPopulation <<
"</UpdateCellPopulation>\n";
742 *rParamsFile <<
"\t\t<SamplingTimestepMultiple>" << mSamplingTimestepMultiple <<
"</SamplingTimestepMultiple>\n";
743 *rParamsFile <<
"\t\t<UpdatingTimestepMultiple>" << mUpdatingTimestepMultiple <<
"</UpdatingTimestepMultiple>\n";
744 *rParamsFile <<
"\t\t<OutputDivisionLocations>" << mOutputDivisionLocations <<
"</OutputDivisionLocations>\n";
745 *rParamsFile <<
"\t\t<OutputCellVelocities>" << mOutputCellVelocities <<
"</OutputCellVelocities>\n";
const double DOUBLE_UNSET
#define EXCEPTION(message)
void SetEndTime(double endTime)
AbstractCellPopulation< ELEMENT_DIM, SPACE_DIM > & mrCellPopulation
void SetUpdateCellPopulationRule(bool updateCellPopulation)
bool GetUpdateCellPopulationRule()
void SetOutputCellVelocities(bool outputCellVelocities)
void SetUpdatingTimestepMultiple(unsigned updatingTimestepMultiple)
virtual bool StoppingEventHasOccurred()
AbstractCellPopulation< ELEMENT_DIM, SPACE_DIM > & rGetCellPopulation()
void RemoveAllCellKillers()
virtual unsigned DoCellBirth()
void SetNoBirth(bool noBirth)
std::vector< boost::shared_ptr< AbstractCellBasedSimulationModifier< ELEMENT_DIM, SPACE_DIM > > > * GetSimulationModifiers()
virtual void UpdateCellPopulation()
std::string GetOutputDirectory()
void AddSimulationModifier(boost::shared_ptr< AbstractCellBasedSimulationModifier< ELEMENT_DIM, SPACE_DIM > > pSimulationModifier)
void AddTopologyUpdateSimulationModifier(boost::shared_ptr< AbstractCellBasedSimulationModifier< ELEMENT_DIM, SPACE_DIM > > pSimulationModifier)
virtual void OutputSimulationParameters(out_stream &rParamsFile)=0
bool GetOutputDivisionLocations()
void OutputSimulationSetup()
bool GetOutputCellVelocities()
void SetOutputDirectory(std::string outputDirectory)
void AddCellKiller(boost::shared_ptr< AbstractCellKiller< SPACE_DIM > > pCellKiller)
virtual ~AbstractCellBasedSimulation()
std::vector< boost::shared_ptr< AbstractCellBasedSimulationModifier< ELEMENT_DIM, SPACE_DIM > > > * GetTopologyUpdateSimulationModifiers()
AbstractCellBasedSimulation(AbstractCellPopulation< ELEMENT_DIM, SPACE_DIM > &rCellPopulation, bool deleteCellPopulationInDestructor=false, bool initialiseCells=true)
std::vector< double > GetNodeLocation(const unsigned &rNodeIndex)
void SetSamplingTimestepMultiple(unsigned samplingTimestepMultiple)
void SetOutputDivisionLocations(bool outputDivisionLocations)
virtual double GetDefaultTimeStep()=0
static void SetOutputDirectory(const std::string &rOutputDirectory)
static void GetBuildInfo(std::string &rInfo)
static void WriteMachineInfoFile(std::string fileBaseName)
static void BeginEvent(unsigned event)
static void EndEvent(unsigned event)
std::string GetOutputDirectoryFullPath() const
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
static RandomNumberGenerator * Instance()
void SetEndTimeAndNumberOfTimeSteps(double endTime, unsigned totalTimeStepsInSimulation)
bool IsEndTimeAndNumberOfTimeStepsSetUp() const
void IncrementTimeOneStep()
static SimulationTime * Instance()
void ResetEndTimeAndNumberOfTimeSteps(const double &rEndTime, const unsigned &rNumberOfTimeStepsInThisRun)
unsigned GetTimeStepsElapsed() const