36#include "PseudoEcgCalculator.hpp"
37#include "HeartRegionCodes.hpp"
38#include "HeartConfig.hpp"
43template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
45 c_vector<double,PROBLEM_DIM>& rU,
46 c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU)
48 c_vector<double,SPACE_DIM> r_vector = rX.
rGetLocation()- mProbeElectrode.rGetLocation();
49 double norm_r = norm_2(r_vector);
50 if (norm_r <= DBL_EPSILON)
52 EXCEPTION(
"Probe is on a mesh Gauss point.");
54 c_vector<double,SPACE_DIM> grad_one_over_r = - (r_vector)*SmallPow( (1/norm_r) , 3);
55 matrix_row<c_matrix<double, PROBLEM_DIM, SPACE_DIM> > grad_u_row(rGradU, 0);
56 double integrand = inner_prod(grad_u_row, grad_one_over_r);
58 return -mDiffusionCoefficient*integrand;
61template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
65 const std::string& rHdf5FileName,
66 const std::string& rVariableName,
67 unsigned timestepStride)
69 mProbeElectrode(rProbeElectrode),
70 mVariableName(rVariableName),
71 mTimestepStride(timestepStride)
81template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
85 if (mVariableName==
"V")
93template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
99template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
102 assert(diffusionCoefficient>=0);
103 mDiffusionCoefficient = diffusionCoefficient;
106template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
110 mpDataReader->GetVariableOverNodes(solution_at_one_time_step, mVariableName , timeStep);
112 double pseudo_ecg_at_one_timestep;
115 pseudo_ecg_at_one_timestep = this->Calculate(mrMesh, solution_at_one_time_step);
123 return pseudo_ecg_at_one_timestep;
126template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
130 std::vector<double> time_points = mpDataReader->GetUnlimitedDimensionValues();
132 std::stringstream stream;
133 stream <<
"PseudoEcgFromElectrodeAt" <<
"_" << mProbeElectrode.GetWithDefault(0)
134 <<
"_" << mProbeElectrode.GetWithDefault(1)
135 <<
"_" << mProbeElectrode.GetWithDefault(2) <<
".dat";
137 out_stream p_file=out_stream(NULL);
141 *p_file <<
"#Time(ms)\tPseudo-ECG\n";
143 for (
unsigned time_step = 0; time_step < mNumTimeSteps; time_step+=mTimestepStride)
145 double pseudo_ecg_at_one_timestep = ComputePseudoEcgAtOneTimeStep(time_step);
148 *p_file << time_points[time_step] <<
"\t" <<pseudo_ecg_at_one_timestep <<
"\n";
#define EXCEPTION(message)
unsigned GetUnsignedAttribute()
static std::string GetProvenanceString()
c_vector< double, DIM > & rGetLocation()
unsigned GetNumberOfRows()
std::vector< double > GetVariableOverTime(const std::string &rVariableName, unsigned nodeIndex)
static HeartConfig * Instance()
static bool IsRegionBath(HeartRegionType regionId)
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
std::string mVariableName
void SetDiffusionCoefficient(double diffusionCoefficient)
bool ShouldSkipThisElement(Element< ELEMENT_DIM, SPACE_DIM > &rElement)
PseudoEcgCalculator(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, const ChastePoint< SPACE_DIM > &rProbeElectrode, const FileFinder &rDirectory, const std::string &rHdf5FileName, const std::string &rVariableName="V", unsigned timestepStride=1)
double ComputePseudoEcgAtOneTimeStep(unsigned timeStep)
Hdf5DataReader * mpDataReader
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > & mrMesh
double mDiffusionCoefficient