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 #ifndef CARDIACSIMULATION_HPP_
00030 #define CARDIACSIMULATION_HPP_
00031
00032
00033 #include "CardiacSimulationArchiver.hpp"
00034
00035
00036
00037
00038
00039
00040 #include <vector>
00041 #include <ctime>
00042 #include <memory>
00043
00044 #include "UblasIncludes.hpp"
00045
00046 #include "MonodomainProblem.hpp"
00047 #include "BidomainProblem.hpp"
00048 #include "BidomainWithBathProblem.hpp"
00049 #include "PetscTools.hpp"
00050 #include "TimeStepper.hpp"
00051 #include "Exception.hpp"
00052
00053 #include "HeartConfig.hpp"
00054 #include "HeartConfigRelatedCellFactory.hpp"
00055 #include "HeartFileFinder.hpp"
00056
00057 #include "TetrahedralMesh.hpp"
00058 #include "NonCachedTetrahedralMesh.hpp"
00059 #include "ChastePoint.hpp"
00060 #include "ChasteCuboid.hpp"
00061 #include "MeshalyzerMeshWriter.hpp"
00062 #include "TrianglesMeshWriter.hpp"
00063
00064 #include "OrthotropicConductivityTensors.hpp"
00065
00066 #include "Hdf5ToMeshalyzerConverter.hpp"
00067 #include "PostProcessingWriter.hpp"
00068
00069 #include "OutputDirectoryFifoQueue.hpp"
00070 #include "ExecutableSupport.hpp"
00071
00082 class CardiacSimulation
00083 {
00084 private:
00090 void ReadParametersFromFile(std::string parameterFileName);
00091
00096 template<class Problem, unsigned SPACE_DIM>
00097 void CreateAndRun()
00098 {
00099 std::auto_ptr<Problem> p_problem;
00100
00101 if (HeartConfig::Instance()->IsSimulationDefined())
00102 {
00103 HeartConfigRelatedCellFactory<SPACE_DIM> cell_factory;
00104 p_problem.reset(new Problem(&cell_factory));
00105
00106 p_problem->Initialise();
00107 }
00108 else
00109 {
00110 p_problem.reset(CardiacSimulationArchiver<Problem>::Load(HeartConfig::Instance()->GetArchivedSimulationDir()));
00111
00112 HeartConfigRelatedCellFactory<SPACE_DIM> cell_factory;
00113 cell_factory.SetMesh(&(p_problem->rGetMesh()));
00114 AbstractCardiacTissue<SPACE_DIM, SPACE_DIM>* p_tissue = p_problem->GetTissue();
00115 DistributedVectorFactory* p_vector_factory = p_problem->rGetMesh().GetDistributedVectorFactory();
00116 for (unsigned node_global_index = p_vector_factory->GetLow();
00117 node_global_index < p_vector_factory->GetHigh();
00118 node_global_index++)
00119 {
00120
00121 cell_factory.SetCellIntracellularStimulus(p_tissue->GetCardiacCell(node_global_index), node_global_index);
00122
00123 cell_factory.SetCellParameters(p_tissue->GetCardiacCell(node_global_index), node_global_index);
00124 }
00125 }
00126
00127 if (HeartConfig::Instance()->GetCheckpointSimulation())
00128 {
00129
00130 OutputDirectoryFifoQueue directory_queue(HeartConfig::Instance()->GetOutputDirectory() + "_checkpoints/",
00131 HeartConfig::Instance()->GetMaxCheckpointsOnDisk());
00132
00133 TimeStepper checkpoint_stepper(p_problem->GetCurrentTime(), HeartConfig::Instance()->GetSimulationDuration(), HeartConfig::Instance()->GetCheckpointTimestep());
00134 while ( !checkpoint_stepper.IsTimeAtEnd() )
00135 {
00136
00137 HeartConfig::Instance()->SetSimulationDuration(checkpoint_stepper.GetNextTime());
00138 p_problem->Solve();
00139
00140
00141 std::stringstream checkpoint_id;
00142 checkpoint_id << HeartConfig::Instance()->GetSimulationDuration() << "ms/";
00143 std::string checkpoint_dir_basename = directory_queue.CreateNextDir(checkpoint_id.str());
00144
00145
00146 std::stringstream archive_foldername;
00147 archive_foldername << HeartConfig::Instance()->GetOutputDirectory() << "_" << HeartConfig::Instance()->GetSimulationDuration() << "ms";
00148 CardiacSimulationArchiver<Problem>::Save(*(p_problem.get()), checkpoint_dir_basename + archive_foldername.str(), false);
00149
00150
00151 OutputFileHandler checkpoint_dir_basename_handler(checkpoint_dir_basename, false);
00152 OutputFileHandler partial_output_dir_handler(HeartConfig::Instance()->GetOutputDirectory(), false);
00153 if (PetscTools::AmMaster())
00154 {
00155 EXPECT0(system, "cp -r " + partial_output_dir_handler.GetOutputDirectoryFullPath() + " " + checkpoint_dir_basename_handler.GetOutputDirectoryFullPath());
00156 }
00157
00158
00159 CreateResumeXmlFile(checkpoint_dir_basename, archive_foldername.str());
00160
00161
00162 checkpoint_stepper.AdvanceOneTimeStep();
00163 }
00164 }
00165 else
00166 {
00167 p_problem->Solve();
00168 }
00169 if (mSaveProblemInstance)
00170 {
00171 mSavedProblem = p_problem;
00172 }
00173 }
00174
00180 void Run();
00181
00191 void CreateResumeXmlFile(const std::string& rOutputDirectory, const std::string& rArchiveDirectory);
00192
00197 std::string BoolToString(bool yesNo);
00198 public:
00208 CardiacSimulation(std::string parameterFileName,
00209 bool writeProvenanceInfo=false,
00210 bool saveProblemInstance=false);
00211
00216 boost::shared_ptr<AbstractUntemplatedCardiacProblem> GetSavedProblem();
00217 private:
00219 bool mSaveProblemInstance;
00220
00222 boost::shared_ptr<AbstractUntemplatedCardiacProblem> mSavedProblem;
00223 };
00224
00225
00226
00227
00228
00229 boost::shared_ptr<AbstractUntemplatedCardiacProblem> CardiacSimulation::GetSavedProblem()
00230 {
00231 return mSavedProblem;
00232 }
00233
00234 std::string CardiacSimulation::BoolToString(bool yesNo)
00235 {
00236 std::string result;
00237 if (yesNo)
00238 {
00239 result = "yes";
00240 }
00241 else
00242 {
00243 result = "no";
00244 }
00245 return result;
00246 }
00247
00248 void CardiacSimulation::CreateResumeXmlFile(const std::string& rOutputDirectory, const std::string& rArchiveDirectory)
00249 {
00250 OutputFileHandler handler(rOutputDirectory, false);
00251 out_stream p_file = handler.OpenOutputFile("ResumeParameters.xml");
00252 (*p_file) << "<?xml version='1.0' encoding='UTF-8'?>" << std::endl;
00253 (*p_file) << "<ChasteParameters xmlns='https://chaste.comlab.ox.ac.uk/nss/parameters/2_1' "
00254 << "xmlns:xsi='http://www.w3.org/2001/XMLSchema-instance' "
00255 << "xsi:schemaLocation='https://chaste.comlab.ox.ac.uk/nss/parameters/2_1 ChasteParameters_2_1.xsd'>" << std::endl;
00256 (*p_file) << std::endl;
00257 (*p_file) << " <ResumeSimulation>" << std::endl;
00258 (*p_file) << " <ArchiveDirectory relative_to='chaste_test_output'>" << rArchiveDirectory << "</ArchiveDirectory>" << std::endl;
00259 (*p_file) << " <SpaceDimension>" << HeartConfig::Instance()->GetSpaceDimension() << "</SpaceDimension>" << std::endl;
00260 (*p_file) << " <SimulationDuration unit='ms'>0.0</SimulationDuration> <!-- Edit with new simulation duration. Please "
00261 << "note that the simulation does not restart at t=0 but at the time where the checkpoint was created.-->" << std::endl;
00262 (*p_file) << " <Domain>" << HeartConfig::Instance()->GetDomain() << "</Domain>" << std::endl;
00263 (*p_file) << " <CheckpointSimulation timestep='" << HeartConfig::Instance()->GetCheckpointTimestep()
00264 << "' unit='ms' max_checkpoints_on_disk='" << HeartConfig::Instance()->GetMaxCheckpointsOnDisk()
00265 << "'/> <!-- This is optional; if not given, the loaded simulation will NOT itself be checkpointed -->" << std::endl;
00266 (*p_file) << " <OutputVisualizer meshalyzer='" << BoolToString(HeartConfig::Instance()->GetVisualizeWithMeshalyzer())
00267 << "' vtk='" << BoolToString(HeartConfig::Instance()->GetVisualizeWithVtk())
00268 << "' cmgui='" << BoolToString(HeartConfig::Instance()->GetVisualizeWithCmgui()) << "'/>" << std::endl;
00269 (*p_file) << " </ResumeSimulation>" << std::endl;
00270 (*p_file) << std::endl;
00271 (*p_file) << " <!-- These elements must exist, but their contents are ignored -->" << std::endl;
00272 (*p_file) << " <Physiological/>" << std::endl;
00273 (*p_file) << " <Numerical/>" << std::endl;
00274 (*p_file) << "</ChasteParameters>" << std::endl;
00275 p_file->close();
00276 HeartConfig::Instance()->CopySchema(handler.GetOutputDirectoryFullPath());
00277 }
00278
00279 CardiacSimulation::CardiacSimulation(std::string parameterFileName,
00280 bool writeProvenanceInfo,
00281 bool saveProblemInstance)
00282 : mSaveProblemInstance(saveProblemInstance)
00283 {
00284
00285 if (parameterFileName == "")
00286 {
00287 EXCEPTION("No XML file name given");
00288 }
00289 ReadParametersFromFile(parameterFileName);
00290 Run();
00291 HeartEventHandler::Headings();
00292 HeartEventHandler::Report();
00293 if (writeProvenanceInfo)
00294 {
00295 ExecutableSupport::SetOutputDirectory(HeartConfig::Instance()->GetOutputDirectory());
00296 ExecutableSupport::WriteProvenanceInfoFile();
00297 ExecutableSupport::WriteMachineInfoFile("machine_info");
00298 }
00299 }
00300
00301 void CardiacSimulation::ReadParametersFromFile(std::string parameterFileName)
00302 {
00303
00304 HeartConfig::Reset();
00305 try
00306 {
00307
00308 HeartConfig::Instance()->SetUseFixedSchemaLocation(true);
00309 HeartConfig::Instance()->SetParametersFile(parameterFileName);
00310 }
00311 catch (Exception& e)
00312 {
00313 if (e.CheckShortMessageContains("Missing file parsing configuration") == "")
00314 {
00315
00316 HeartConfig::Reset();
00317 HeartConfig::Instance()->SetUseFixedSchemaLocation(false);
00318 HeartConfig::Instance()->SetParametersFile(parameterFileName);
00319 }
00320 else
00321 {
00322 throw e;
00323 }
00324 }
00325 }
00326
00327
00328 #define DOMAIN_CASE(VALUE, CLASS, DIM) \
00329 case VALUE: \
00330 { \
00331 CreateAndRun<CLASS<DIM>, DIM>(); \
00332 break; \
00333 }
00334
00335 #define DOMAIN_SWITCH(DIM) \
00336 switch (HeartConfig::Instance()->GetDomain()) \
00337 { \
00338 DOMAIN_CASE(cp::domain_type::Mono, MonodomainProblem, DIM) \
00339 DOMAIN_CASE(cp::domain_type::Bi, BidomainProblem, DIM) \
00340 DOMAIN_CASE(cp::domain_type::BiWithBath, BidomainWithBathProblem, DIM) \
00341 default: \
00342 NEVER_REACHED; \
00343 } \
00344 break
00345
00346
00347
00348 void CardiacSimulation::Run()
00349 {
00350 switch (HeartConfig::Instance()->GetSpaceDimension())
00351 {
00352 case 3:
00353 {
00354 DOMAIN_SWITCH(3);
00355 }
00356 case 2:
00357 {
00358 DOMAIN_SWITCH(2);
00359 }
00360 case 1:
00361 {
00362 DOMAIN_SWITCH(1);
00363 }
00364 default:
00365
00366 EXCEPTION("Space dimension not supported: should be 1, 2 or 3");
00367 }
00368 }
00369
00370
00371 #undef DOMAIN_SWITCH
00372 #undef DOMAIN_CASE
00373
00374
00375 #endif