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 #include "PseudoEcgCalculator.hpp"
00030 #include "HeartConfig.hpp"
00031 #include "PetscTools.hpp"
00032 #include "Version.hpp"
00033 #include <iostream>
00034
00035 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00036 double PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> ::GetIntegrand(ChastePoint<SPACE_DIM>& rX,
00037 c_vector<double,PROBLEM_DIM>& rU,
00038 c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU)
00039 {
00040 c_vector<double,SPACE_DIM> r_vector = rX.rGetLocation()- mProbeElectrode.rGetLocation();
00041 c_vector<double,SPACE_DIM> grad_one_over_r = - (r_vector)*SmallPow( (1/norm_2(r_vector)) , 3);
00042 matrix_row<c_matrix<double, PROBLEM_DIM, SPACE_DIM> > grad_u_row(rGradU, 0);
00043 double integrand = inner_prod(grad_u_row, grad_one_over_r);
00044
00045 return -mDiffusionCoefficient*integrand;
00046 }
00047
00048 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00049 PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> ::PseudoEcgCalculator (AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>& rMesh,
00050 const ChastePoint<SPACE_DIM>& rProbeElectrode,
00051 std::string directory,
00052 std::string hdf5File,
00053 std::string variableName,
00054 bool makeAbsolute)
00055 : mrMesh(rMesh),
00056 mProbeElectrode(rProbeElectrode),
00057 mVariableName(variableName)
00058
00059 {
00060 mpDataReader = new Hdf5DataReader(directory, hdf5File, makeAbsolute);
00061 mNumberOfNodes = mpDataReader->GetNumberOfRows();
00062 mNumTimeSteps = mpDataReader->GetVariableOverTime(mVariableName, 0u).size();
00063 mDiffusionCoefficient = 1.0;
00064
00065 assert(mNumberOfNodes == mrMesh.GetNumNodes());
00066 }
00067
00068 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00069 PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::~PseudoEcgCalculator()
00070 {
00071 delete mpDataReader;
00072 }
00073
00074 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00075 void PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetDiffusionCoefficient(double diffusionCoefficient)
00076 {
00077 assert(diffusionCoefficient>=0);
00078 mDiffusionCoefficient = diffusionCoefficient;
00079 }
00080
00081 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00082 double PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ComputePseudoEcgAtOneTimeStep (unsigned timeStep)
00083 {
00084 Vec solution_at_one_time_step = PetscTools::CreateVec(mNumberOfNodes);
00085 mpDataReader->GetVariableOverNodes(solution_at_one_time_step, mVariableName , timeStep);
00086
00087 double pseudo_ecg_at_one_timestep = Calculate(mrMesh, solution_at_one_time_step);
00088
00089 VecDestroy(solution_at_one_time_step);
00090 return pseudo_ecg_at_one_timestep;
00091 }
00092
00093 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00094 void PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::WritePseudoEcg ()
00095 {
00096 std::stringstream stream;
00097 stream << "PseudoEcg.dat";
00098 OutputFileHandler output_file_handler(HeartConfig::Instance()->GetOutputDirectory() + "/output", false);
00099 out_stream p_file=out_stream(NULL);
00100 if (PetscTools::AmMaster())
00101 {
00102 p_file = output_file_handler.OpenOutputFile(stream.str());
00103 }
00104 for (unsigned i = 0; i < mNumTimeSteps; i++)
00105 {
00106 double pseudo_ecg_at_one_timestep = ComputePseudoEcgAtOneTimeStep(i);
00107 if (PetscTools::AmMaster())
00108 {
00109 *p_file << pseudo_ecg_at_one_timestep << "\n";
00110 }
00111 }
00112 if (PetscTools::AmMaster())
00113 {
00114
00115 std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
00116 *p_file << comment;
00117 p_file->close();
00118 }
00119 }
00121
00123
00124 template class PseudoEcgCalculator<1,1,1>;
00125
00126
00127 template class PseudoEcgCalculator<2,2,1>;
00128
00129 template class PseudoEcgCalculator<3,3,1>;
00130