35#include "AbstractCardiacProblem.hpp"
37#include "DistributedVector.hpp"
39#include "GenericMeshReader.hpp"
40#include "Hdf5ToCmguiConverter.hpp"
41#include "Hdf5ToMeshalyzerConverter.hpp"
42#include "Hdf5ToVtkConverter.hpp"
43#include "HeartConfig.hpp"
44#include "HeartEventHandler.hpp"
45#include "LinearSystem.hpp"
47#include "PostProcessingWriter.hpp"
48#include "ProgressReporter.hpp"
49#include "TimeStepper.hpp"
51template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
55 mAllocatedMemoryForMesh(false),
58 mpCardiacTissue(NULL),
60 mpCellFactory(pCellFactory),
64 mpTimeAdaptivityController(NULL),
66 mUseHdf5DataWriterCache(false),
67 mHdf5DataWriterChunkSizeAndAlignment(0)
72 EXCEPTION(
"AbstractCardiacProblem: Please supply a cell factory pointer to your cardiac problem constructor.");
77template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
82 mAllocatedMemoryForMesh(false),
85 mVoltageColumnId(UINT_MAX),
86 mTimeColumnId(UINT_MAX),
87 mNodeColumnId(UINT_MAX),
88 mpCardiacTissue(NULL),
94 mpTimeAdaptivityController(NULL),
96 mUseHdf5DataWriterCache(false),
97 mHdf5DataWriterChunkSizeAndAlignment(0)
101template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
104 delete mpCardiacTissue;
110 if (mAllocatedMemoryForMesh)
116template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
124 WARNING(
"Using a non-distributed mesh in a parallel simulation is not a good idea.");
134 CreateMeshFromHeartConfig();
135 std::shared_ptr<AbstractMeshReader<ELEMENT_DIM, SPACE_DIM> > p_mesh_reader
137 mpMesh->ConstructFromMeshReader(*p_mesh_reader);
141 CreateMeshFromHeartConfig();
149 c_vector<double, 1> fibre_length;
151 mpMesh->ConstructRegularSlabMesh(inter_node_space, fibre_length[0]);
156 c_vector<double, 2> sheet_dimensions;
158 mpMesh->ConstructRegularSlabMesh(inter_node_space, sheet_dimensions[0], sheet_dimensions[1]);
163 c_vector<double, 3> slab_dimensions;
165 mpMesh->ConstructRegularSlabMesh(inter_node_space, slab_dimensions[0], slab_dimensions[1], slab_dimensions[2]);
177 mAllocatedMemoryForMesh =
true;
181 EXCEPTION(std::string(
"No mesh given: define it in XML parameters file or call SetMesh()\n") + e.GetShortMessage());
184 mpCellFactory->SetMesh(mpMesh);
192 mpCellFactory->FillInCellularTransmuralAreas();
195 delete mpCardiacTissue;
196 mpCardiacTissue = CreateCardiacTissue();
216template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
222template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
225 this->mpBoundaryConditionsContainer = pBcc;
228template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
231 if (mpCardiacTissue == NULL)
233 EXCEPTION(
"Cardiac tissue is null, Initialise() probably hasn't been called");
237 EXCEPTION(
"End time should be in the future");
243 EXCEPTION(
"Either explicitly specify not to print output (call PrintOutput(false)) or specify the output directory and filename prefix");
258 if (fabs(end_time - pde_time * round(end_time / pde_time)) > 1e-10)
260 EXCEPTION(
"PDE timestep does not seem to divide end time - check parameters");
264template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
268 Vec initial_condition = p_factory->
CreateVec(PROBLEM_DIM);
270 std::vector<DistributedVector::Stripe> stripe;
271 stripe.reserve(PROBLEM_DIM);
273 for (
unsigned i = 0; i < PROBLEM_DIM; i++)
282 stripe[0][index] = mpCardiacTissue->GetCardiacCell(index.Global)->GetVoltage();
283 if (PROBLEM_DIM == 2)
285 stripe[1][index] = 0;
291 return initial_condition;
294template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
301 assert(mpMesh == NULL);
302 assert(pMesh != NULL);
303 mAllocatedMemoryForMesh =
false;
307template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
310 mPrintOutput = printOutput;
313template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
316 mWriteInfo = writeInfo;
319template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
325template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
328 return mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(mSolution);
331template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
337template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
344template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
347 if (mpCardiacTissue == NULL)
349 EXCEPTION(
"Tissue not yet set up, you may need to call Initialise() before GetTissue().");
351 return mpCardiacTissue;
354template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
362 mpTimeAdaptivityController = pController;
366 mpTimeAdaptivityController = NULL;
370template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
375 std::vector<double> additional_stopping_times;
376 SetUpAdditionalStoppingTimes(additional_stopping_times);
382 additional_stopping_times);
387 if (!mpBoundaryConditionsContainer)
391 for (
unsigned problem_index = 0; problem_index < PROBLEM_DIM; problem_index++)
393 mpDefaultBoundaryConditionsContainer->DefineZeroNeumannOnMeshBoundary(mpMesh, problem_index);
395 mpBoundaryConditionsContainer = mpDefaultBoundaryConditionsContainer;
398 assert(mpSolver == NULL);
399 mpSolver = CreateSolver();
402 Vec initial_condition;
405 initial_condition = mSolution;
409 initial_condition = CreateInitialCondition();
412 std::string progress_reporter_dir;
417 bool extending_file =
false;
420 extending_file = InitialiseWriter();
427 if (mSolution != initial_condition)
447 if (!(mSolution && extending_file))
449 WriteOneStep(stepper.
GetTime(), initial_condition);
450 mpWriter->AdvanceAlongUnlimitedDimension();
458 progress_reporter_dir =
"";
460 BOOST_FOREACH (boost::shared_ptr<AbstractOutputModifier> p_output_modifier, mOutputModifiers)
462 p_output_modifier->InitialiseAtStart(this->mpMesh->GetDistributedVectorFactory(), this->mpMesh->rGetNodePermutation());
463 p_output_modifier->ProcessSolutionAtTimeStep(stepper.
GetTime(), initial_condition, PROBLEM_DIM);
474 progress_reporter.Update(mCurrentTime);
477 if (mpTimeAdaptivityController)
479 mpSolver->SetTimeAdaptivityController(mpTimeAdaptivityController);
486 mpSolver->SetInitialCondition(initial_condition);
488 AtBeginningOfTimestep(stepper.
GetTime());
494 mSolution = mpSolver->Solve();
512 if (initial_condition != mSolution)
527 CloseFilesAndPostProcess();
538 initial_condition = mSolution;
542 mCurrentTime = stepper.
GetTime();
552 BOOST_FOREACH (boost::shared_ptr<AbstractOutputModifier> p_output_modifier, mOutputModifiers)
554 p_output_modifier->ProcessSolutionAtTimeStep(stepper.
GetTime(), mSolution, PROBLEM_DIM);
560 WriteOneStep(stepper.
GetTime(), mSolution);
562 mpWriter->AdvanceAlongUnlimitedDimension();
567 progress_reporter.Update(stepper.
GetTime());
569 OnEndOfTimestep(stepper.
GetTime());
577 progress_reporter.PrintFinalising();
578 BOOST_FOREACH (boost::shared_ptr<AbstractOutputModifier> p_output_modifier, mOutputModifiers)
580 p_output_modifier->FinaliseAtEnd();
582 CloseFilesAndPostProcess();
586template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
619 mHdf5DataWriterChunkSizeAndAlignment);
620 post_writer.WritePostProcessingFiles();
630 if (mNodesToOutput.empty())
640 std::string subdirectory_name = converter.GetSubdirectory();
652 std::string subdirectory_name = converter.GetSubdirectory();
664 std::string subdirectory_name = converter.GetSubdirectory();
676 std::string subdirectory_name = converter.GetSubdirectory();
683template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
688 if (mNodesToOutput.empty())
691 mpWriter->DefineFixedDimension(mpMesh->GetNumNodes());
696 if (mpMesh->rGetNodePermutation().size() > 0)
700 EXCEPTION(
"HeartConfig setting `GetOutputUsingOriginalNodeOrdering` is meaningless when outputting particular nodes in parallel. (Nodes are written with their original indices by default).");
702 std::vector<unsigned> nodes_to_output_permuted(mNodesToOutput.size());
703 for (
unsigned i = 0; i < mNodesToOutput.size(); i++)
705 nodes_to_output_permuted[i] = mpMesh->rGetNodePermutation()[mNodesToOutput[i]];
707 mpWriter->DefineFixedDimension(mNodesToOutput, nodes_to_output_permuted, mpMesh->GetNumNodes());
710 mpWriter->DefineFixedDimension(mNodesToOutput, mNodesToOutput, mpMesh->GetNumNodes());
714 mVoltageColumnId = mpWriter->DefineVariable(
"V",
"mV");
721 mpWriter->DefineUnlimitedDimension(
"Time",
"msecs", stepper.
EstimateTimeSteps() + 1);
725 mVoltageColumnId = mpWriter->GetVariableByName(
"V");
729template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
732 mExtraVariablesId.clear();
737 std::vector<std::string> output_variables;
739 const unsigned num_vars = output_variables.size();
740 mExtraVariablesId.reserve(num_vars);
743 for (
unsigned var_index = 0; var_index < num_vars; var_index++)
746 std::string var_name = output_variables[var_index];
752 column_id = this->mpWriter->GetVariableByName(var_name);
758 column_id = this->mpWriter->DefineVariable(var_name,
"unknown_units");
762 mExtraVariablesId.push_back(column_id);
767template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
771 std::vector<std::string> output_variables;
772 unsigned num_vars = mExtraVariablesId.size();
777 assert(output_variables.size() == num_vars);
780 for (
unsigned var_index = 0; var_index < num_vars; var_index++)
783 Vec variable_data = this->mpMesh->GetDistributedVectorFactory()->CreateVec();
784 DistributedVector distributed_var_data = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(variable_data);
788 index != distributed_var_data.
End();
796 distributed_var_data[index] = 0.0;
801 distributed_var_data[index] = this->mpCardiacTissue->GetCardiacCell(index.Global)->GetAnyVariable(output_variables[var_index], mCurrentTime);
804 distributed_var_data.
Restore();
807 this->mpWriter->PutVector(mExtraVariablesId[var_index], variable_data);
813template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
816 bool extend_file = (mSolution != NULL);
841 if (times.back() > mCurrentTime)
843 EXCEPTION(
"Attempting to extend " << h5_file.
GetAbsolutePath() <<
" with results from time = " << mCurrentTime <<
", but it already contains results up to time = " << times.back() <<
"."
844 " Calling HeartConfig::Instance()->SetOutputDirectory() before Solve() will direct results elsewhere.");
849 mpWriter =
new Hdf5DataWriter(*mpMesh->GetDistributedVectorFactory(),
855 mUseHdf5DataWriterCache);
864 if (!extend_file && mHdf5DataWriterChunkSizeAndAlignment)
866 mpWriter->SetTargetChunkSize(mHdf5DataWriterChunkSizeAndAlignment);
867 mpWriter->SetAlignment(mHdf5DataWriterChunkSizeAndAlignment);
871 DefineWriterColumns(extend_file);
876 bool success = mpWriter->ApplyPermutation(mpMesh->rGetNodePermutation(),
true );
877 if (success ==
false)
886 mpWriter->EndDefineMode();
892template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
895 mUseHdf5DataWriterCache = useCache;
898template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
901 mHdf5DataWriterChunkSizeAndAlignment = size;
904template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
907 mNodesToOutput = nodesToOutput;
910template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
915 EXCEPTION(
"Data reader invalid as data writer cannot be initialised");
920template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
926template <
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
#define EXCEPTION(message)
virtual Vec CreateInitialCondition()
virtual ~AbstractCardiacProblem()
void SetWriteInfo(bool writeInfo=true)
virtual void SetElectrodes()
virtual void CreateMeshFromHeartConfig()
void SetOutputNodes(std::vector< unsigned > &rNodesToOutput)
void DefineExtraVariablesWriterColumns(bool extending)
DistributedVector GetSolutionDistributedVector()
void SetHdf5DataWriterTargetChunkSizeAndAlignment(hsize_t size)
Hdf5DataReader GetDataReader()
void SetUseTimeAdaptivityController(bool useAdaptivity, AbstractTimeAdaptivityController *pController=NULL)
void SetMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > & rGetMesh()
void CloseFilesAndPostProcess()
virtual void PreSolveChecks()
void PrintOutput(bool rPrintOutput)
virtual bool GetHasBath()
void SetBoundaryConditionsContainer(BccType pBcc)
AbstractCardiacCellFactory< ELEMENT_DIM, SPACE_DIM > * mpCellFactory
std::vector< unsigned > mNodesToOutput
AbstractCardiacTissue< ELEMENT_DIM, SPACE_DIM > * GetTissue()
void SetUseHdf5DataWriterCache(bool useCache=true)
void WriteExtraVariablesOneStep()
virtual void DefineWriterColumns(bool extending)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
std::string GetAbsolutePath() const
static void BeginEvent(unsigned event)
static void EndEvent(unsigned event)
std::vector< double > GetUnlimitedDimensionValues()
void GetSheetDimensions(c_vector< double, 2 > &sheetDimensions) const
double GetPdeTimeStep() const
void GetOutputVariables(std::vector< std::string > &rOutputVariables) const
double GetSimulationDuration() const
void SetOutputUsingOriginalNodeOrdering(bool useOriginal)
std::string GetOutputFilenamePrefix() const
void GetSlabDimensions(c_vector< double, 3 > &slabDimensions) const
void Write(bool useArchiveLocationInfo=false, std::string subfolderName="output")
double GetInterNodeSpace() const
std::string GetOutputDirectory() const
void GetFibreLength(c_vector< double, 1 > &fibreLength) const
static HeartConfig * Instance()
static bool IsRegionBath(HeartRegionType regionId)
static std::string GetChasteTestOutputDirectory()
void AdvanceOneTimeStep()
double GetNextTime() const
unsigned EstimateTimeSteps() const