37#include "HeartConfig.hpp"
38#include "PostProcessingWriter.hpp"
40#include "OutputFileHandler.hpp"
41#include "DistanceMapCalculator.hpp"
42#include "PseudoEcgCalculator.hpp"
44#include "HeartEventHandler.hpp"
45#include "Hdf5DataWriter.hpp"
46#include "Hdf5ToMeshalyzerConverter.hpp"
47#include "Hdf5ToVtkConverter.hpp"
51template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
54 const std::string& rHdf5FileName,
55 const std::string& rVoltageName,
56 hsize_t hdf5DataWriterChunkSize)
57 : mDirectory(rDirectory),
58 mHdf5File(rHdf5FileName),
59 mVoltageName(rVoltageName),
61 mHdf5DataWriterChunkSize(hdf5DataWriterChunkSize)
63 mLo =
mrMesh.GetDistributedVectorFactory()->GetLow();
64 mHi =
mrMesh.GetDistributedVectorFactory()->GetHigh();
71template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
81 std::vector<std::pair<double,double> > apd_maps;
83 for (
unsigned i=0; i<apd_maps.size(); i++)
85 WriteApdMapFile(apd_maps[i].first, apd_maps[i].second);
91 std::vector<double> upstroke_time_maps;
93 for (
unsigned i=0; i<upstroke_time_maps.size(); i++)
95 WriteUpstrokeTimeMap(upstroke_time_maps[i]);
101 std::vector<double> upstroke_velocity_maps;
103 for (
unsigned i=0; i<upstroke_velocity_maps.size(); i++)
105 WriteMaxUpstrokeVelocityMap(upstroke_velocity_maps[i]);
111 std::vector<unsigned> conduction_velocity_maps;
117 for (
unsigned i=0; i<conduction_velocity_maps.size(); i++)
119 std::vector<double> distance_map;
120 std::vector<unsigned> origin_surface;
121 origin_surface.push_back(conduction_velocity_maps[i]);
123 WriteConductionVelocityMap(conduction_velocity_maps[i], distance_map);
129 std::vector<unsigned> requested_nodes;
131 WriteVariablesOverTimeAtNodes(requested_nodes);
136 std::vector<ChastePoint<SPACE_DIM> > electrodes;
143 for (
unsigned i=0; i<electrodes.size(); i++)
155 mpCalculator->SetHdf5DataReader(mpDataReader);
159template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
166template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
168 const std::string& rDatasetName,
169 const std::string& rDatasetUnit,
170 const std::string& rUnlimitedVariableName,
171 const std::string& rUnlimitedVariableUnit)
176 mDirectory.GetRelativePath(test_output),
191 if (mHdf5DataWriterChunkSize>0u)
208 unsigned local_max_paces = 0u;
209 for (
unsigned node_index = 0; node_index < rDataPayload.size(); ++node_index)
211 if (rDataPayload[node_index].size() > local_max_paces)
213 local_max_paces = rDataPayload[node_index].size();
217 unsigned max_paces = 0u;
218 MPI_Allreduce(&local_max_paces, &max_paces, 1, MPI_UNSIGNED, MPI_MAX, PETSC_COMM_WORLD);
220 for (
unsigned pace_idx = 0; pace_idx < max_paces; pace_idx++)
225 index!= distributed_vector.
End();
228 unsigned node_idx = index.Local;
230 if (pace_idx < rDataPayload[node_idx].size() )
232 distributed_vector[index] = rDataPayload[node_idx][pace_idx];
236 distributed_vector[index] = -999.0;
247template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
250 std::vector<std::vector<double> > local_output_data = mpCalculator->CalculateAllActionPotentialDurationsForNodeRange(repolarisationPercentage, mLo, mHi, threshold);
253 std::stringstream hdf5_dataset_name;
254 hdf5_dataset_name <<
"Apd_" << repolarisationPercentage;
256 WriteOutputDataToHdf5(local_output_data,
257 hdf5_dataset_name.str() + ConvertToHdf5FriendlyString(threshold) +
"_Map",
261template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
264 std::vector<std::vector<double> > output_data;
266 for (
unsigned node_index = mLo; node_index < mHi; node_index++)
268 std::vector<double> upstroke_times;
271 upstroke_times = mpCalculator->CalculateUpstrokeTimes(node_index, threshold);
272 assert(upstroke_times.size() != 0);
276 upstroke_times.push_back(0);
277 assert(upstroke_times.size() == 1);
279 output_data.push_back(upstroke_times);
282 WriteOutputDataToHdf5(output_data,
283 "UpstrokeTimeMap" + ConvertToHdf5FriendlyString(threshold),
287template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
290 std::vector<std::vector<double> > output_data;
292 for (
unsigned node_index = mLo; node_index < mHi; node_index++)
294 std::vector<double> upstroke_velocities;
297 upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
298 assert(upstroke_velocities.size() != 0);
302 upstroke_velocities.push_back(0);
303 assert(upstroke_velocities.size() ==1);
305 output_data.push_back(upstroke_velocities);
308 WriteOutputDataToHdf5(output_data,
309 "MaxUpstrokeVelocityMap" + ConvertToHdf5FriendlyString(threshold),
313template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
317 std::vector<std::vector<double> > output_data;
318 for (
unsigned dest_node = mLo; dest_node < mHi; dest_node++)
320 std::vector<double> conduction_velocities;
323 conduction_velocities = mpCalculator->CalculateAllConductionVelocities(originNode, dest_node, distancesFromOriginNode[dest_node]);
324 assert(conduction_velocities.size() != 0);
328 conduction_velocities.push_back(0);
329 assert(conduction_velocities.size() == 1);
331 output_data.push_back(conduction_velocities);
333 std::stringstream filename_stream;
334 filename_stream <<
"ConductionVelocityFromNode" << originNode;
336 WriteOutputDataToHdf5(output_data, filename_stream.str(),
"cm_per_msec");
339template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
342 std::vector<std::vector<double> > output_data;
345 for (
unsigned node_index = mLo; node_index < mHi; node_index++)
347 std::vector<double> upstroke_velocities;
348 std::vector<unsigned> above_threshold_depolarisations;
349 std::vector<double> output_item;
350 bool no_upstroke_occurred =
false;
354 upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
355 assert(upstroke_velocities.size() != 0);
359 upstroke_velocities.push_back(0);
360 assert(upstroke_velocities.size() == 1);
361 no_upstroke_occurred =
true;
364 above_threshold_depolarisations = mpCalculator->CalculateAllAboveThresholdDepolarisations(node_index, threshold);
367 unsigned total_number_of_above_threshold_depolarisations = 0;
368 for (
unsigned ead_index = 0; ead_index< above_threshold_depolarisations.size();ead_index++)
370 total_number_of_above_threshold_depolarisations = total_number_of_above_threshold_depolarisations + above_threshold_depolarisations[ead_index];
374 if (no_upstroke_occurred)
376 output_item.push_back(0);
380 output_item.push_back((
double)upstroke_velocities.size());
383 output_item.push_back((
double) total_number_of_above_threshold_depolarisations);
385 output_data.push_back(output_item);
389 WriteGenericFileToMeshalyzer(output_data,
"",
"AboveThresholdDepolarisations" + ConvertToHdf5FriendlyString(threshold) +
393template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
396 std::stringstream dataset_stream;
397 std::string extra_message =
"_";
400 extra_message +=
"minus_";
401 threshold = -threshold;
405 if (threshold - floor(threshold) > 1e-8)
408 dataset_stream << extra_message << floor(threshold) <<
"pt" << floor(threshold*100)-(floor(threshold)*100);
412 dataset_stream << extra_message << threshold;
415 return dataset_stream.str();
418template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
421 std::vector<std::string> variable_names = mpDataReader->GetVariableNames();
424 for (
unsigned name_index=0; name_index < variable_names.size(); name_index++)
426 std::vector<std::vector<double> > output_data;
430 output_data.resize( mpDataReader->GetUnlimitedDimensionValues().size() );
431 for (
unsigned j = 0; j < mpDataReader->GetUnlimitedDimensionValues().size(); j++)
433 output_data[j].resize(rNodeIndices.size());
436 for (
unsigned requested_index = 0; requested_index < rNodeIndices.size(); requested_index++)
438 unsigned node_index = rNodeIndices[requested_index];
441 if ((mrMesh.rGetNodePermutation().size() != 0) &&
444 node_index = mrMesh.rGetNodePermutation()[ rNodeIndices[requested_index] ];
448 std::vector<double> time_series = mpDataReader->GetVariableOverTime(variable_names[name_index], node_index);
449 assert ( time_series.size() == mpDataReader->GetUnlimitedDimensionValues().size());
452 for (
unsigned time_step = 0; time_step < time_series.size(); time_step++)
454 output_data[time_step][requested_index] = time_series[time_step];
458 std::stringstream filename_stream;
459 filename_stream <<
"NodalTraces_" << variable_names[name_index] <<
".dat";
461 WriteGenericFileToMeshalyzer(output_data,
"", filename_stream.str());
465template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM>
471 out_stream p_file=out_stream(NULL);
481 p_file = output_file_handler.
OpenOutputFile(rFileName, std::ios::app);
484 for (
unsigned line_number=0; line_number<rDataPayload.size(); line_number++)
486 for (
unsigned i = 0; i < rDataPayload[line_number].size(); i++)
488 *p_file << rDataPayload[line_number][i] <<
"\t";
490 *p_file << std::endl;
static std::string GetProvenanceString()
void ComputeDistanceMap(const std::vector< unsigned > &rSourceNodeIndices, std::vector< double > &rNodeDistances)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
unsigned GetNumberOfRows()
void SetTargetChunkSize(hsize_t targetSize)
void DefineUnlimitedDimension(const std::string &rVariableName, const std::string &rVariableUnits, unsigned estimatedLength=1)
void PutUnlimitedVariable(double value)
void AdvanceAlongUnlimitedDimension()
int GetVariableByName(const std::string &rVariableName)
void DefineFixedDimension(long dimensionSize)
int DefineVariable(const std::string &rVariableName, const std::string &rVariableUnits)
virtual void EndDefineMode()
void PutVector(int variableID, Vec petscVector)
bool GetOutputUsingOriginalNodeOrdering()
void GetNodalTimeTraceRequested(std::vector< unsigned > &rRequestedNodes) const
void GetApdMaps(std::vector< std::pair< double, double > > &rApdMaps) const
void GetConductionVelocityMaps(std::vector< unsigned > &rConductionVelocityMaps) const
void GetPseudoEcgElectrodePositions(std::vector< ChastePoint< SPACE_DIM > > &rPseudoEcgElectrodePositions) const
void GetUpstrokeTimeMaps(std::vector< double > &rUpstrokeTimeMaps) const
void GetMaxUpstrokeVelocityMaps(std::vector< double > &rUpstrokeVelocityMaps) const
static HeartConfig * Instance()
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
void WritePostProcessingFiles()
void WriteMaxUpstrokeVelocityMap(double threshold)
void WriteConductionVelocityMap(unsigned originNode, std::vector< double > distancesFromOriginNode)
void WriteVariablesOverTimeAtNodes(std::vector< unsigned > &rNodeIndices)
void WriteUpstrokeTimeMap(double threshold)
void WriteGenericFileToMeshalyzer(std::vector< std::vector< double > > &rDataPayload, const std::string &rFolder, const std::string &rFileName)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > & mrMesh
std::string ConvertToHdf5FriendlyString(double threshold)
void WriteApdMapFile(double repolarisationPercentage, double threshold)
void WriteOutputDataToHdf5(const std::vector< std::vector< double > > &rDataPayload, const std::string &rDatasetName, const std::string &rDatasetUnit, const std::string &rUnlimitedVariableName="PaceNumber", const std::string &rUnlimitedVariableUnit="dimensionless")
Hdf5DataReader * mpDataReader
PostProcessingWriter(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, const FileFinder &rDirectory, const std::string &rHdf5FileName, const std::string &rVoltageName="V", hsize_t hdf5DataWriterChunkSize=0)
void WriteAboveThresholdDepolarisationFile(double threshold)
PropagationPropertiesCalculator * mpCalculator