Chaste Release::3.1
|
00001 /* 00002 00003 Copyright (c) 2005-2012, University of Oxford. 00004 All rights reserved. 00005 00006 University of Oxford means the Chancellor, Masters and Scholars of the 00007 University of Oxford, having an administrative office at Wellington 00008 Square, Oxford OX1 2JD, UK. 00009 00010 This file is part of Chaste. 00011 00012 Redistribution and use in source and binary forms, with or without 00013 modification, are permitted provided that the following conditions are met: 00014 * Redistributions of source code must retain the above copyright notice, 00015 this list of conditions and the following disclaimer. 00016 * Redistributions in binary form must reproduce the above copyright notice, 00017 this list of conditions and the following disclaimer in the documentation 00018 and/or other materials provided with the distribution. 00019 * Neither the name of the University of Oxford nor the names of its 00020 contributors may be used to endorse or promote products derived from this 00021 software without specific prior written permission. 00022 00023 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 00024 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 00025 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 00026 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE 00027 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 00028 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE 00029 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) 00030 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 00031 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 00032 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00033 00034 */ 00035 00036 #include "PseudoEcgCalculator.hpp" 00037 #include "HeartRegionCodes.hpp" 00038 #include "HeartConfig.hpp" 00039 #include "PetscTools.hpp" 00040 #include "Version.hpp" 00041 #include <iostream> 00042 00043 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM> 00044 double PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> ::GetIntegrand(ChastePoint<SPACE_DIM>& rX, 00045 c_vector<double,PROBLEM_DIM>& rU, 00046 c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU) 00047 { 00048 c_vector<double,SPACE_DIM> r_vector = rX.rGetLocation()- mProbeElectrode.rGetLocation(); 00049 double norm_r = norm_2(r_vector); 00050 if (norm_r <= DBL_EPSILON) 00051 { 00052 EXCEPTION("Probe is on a mesh Gauss point."); 00053 } 00054 c_vector<double,SPACE_DIM> grad_one_over_r = - (r_vector)*SmallPow( (1/norm_r) , 3); 00055 matrix_row<c_matrix<double, PROBLEM_DIM, SPACE_DIM> > grad_u_row(rGradU, 0); 00056 double integrand = inner_prod(grad_u_row, grad_one_over_r); 00057 00058 return -mDiffusionCoefficient*integrand; 00059 } 00060 00061 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM> 00062 PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::PseudoEcgCalculator(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>& rMesh, 00063 const ChastePoint<SPACE_DIM>& rProbeElectrode, 00064 std::string directory, 00065 std::string hdf5File, 00066 std::string variableName, 00067 bool makeAbsolute) 00068 : mrMesh(rMesh), 00069 mProbeElectrode(rProbeElectrode), 00070 mVariableName(variableName) 00071 { 00072 mpDataReader = new Hdf5DataReader(directory, hdf5File, makeAbsolute); 00073 mNumberOfNodes = mpDataReader->GetNumberOfRows(); 00074 mNumTimeSteps = mpDataReader->GetVariableOverTime(mVariableName, 0u).size(); 00075 mDiffusionCoefficient = 1.0; 00076 //check that the hdf file was generated by simulations from the same mesh 00077 assert(mNumberOfNodes == mrMesh.GetNumNodes()); 00078 } 00079 00080 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM> 00081 bool PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ShouldSkipThisElement(Element<ELEMENT_DIM,SPACE_DIM>& rElement) 00082 { 00083 bool result = false; // Don't skip usually 00084 if (mVariableName=="V") 00085 { 00086 // If we are integrating voltage and this element is in the bath then skip it. 00087 result = HeartRegionCode::IsRegionBath(rElement.GetUnsignedAttribute()); 00088 } 00089 return result; 00090 } 00091 00092 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM> 00093 PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::~PseudoEcgCalculator() 00094 { 00095 delete mpDataReader; 00096 } 00097 00098 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM> 00099 void PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetDiffusionCoefficient(double diffusionCoefficient) 00100 { 00101 assert(diffusionCoefficient>=0); 00102 mDiffusionCoefficient = diffusionCoefficient; 00103 } 00104 00105 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM> 00106 double PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::ComputePseudoEcgAtOneTimeStep (unsigned timeStep) 00107 { 00108 Vec solution_at_one_time_step = PetscTools::CreateVec(mNumberOfNodes); 00109 mpDataReader->GetVariableOverNodes(solution_at_one_time_step, mVariableName , timeStep); 00110 00111 double pseudo_ecg_at_one_timestep; 00112 try 00113 { 00114 pseudo_ecg_at_one_timestep = Calculate(mrMesh, solution_at_one_time_step); 00115 } 00116 catch (Exception &e) 00117 { 00118 PetscTools::Destroy(solution_at_one_time_step); 00119 throw e; 00120 } 00121 PetscTools::Destroy(solution_at_one_time_step); 00122 return pseudo_ecg_at_one_timestep; 00123 } 00124 00125 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM> 00126 void PseudoEcgCalculator<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::WritePseudoEcg () 00127 { 00128 std::stringstream stream; 00129 stream << "PseudoEcgFromElectrodeAt" << "_" << mProbeElectrode.GetWithDefault(0) 00130 << "_" << mProbeElectrode.GetWithDefault(1) 00131 << "_" << mProbeElectrode.GetWithDefault(2) << ".dat"; 00132 OutputFileHandler output_file_handler(HeartConfig::Instance()->GetOutputDirectory() + "/output", false); 00133 out_stream p_file=out_stream(NULL); 00134 if (PetscTools::AmMaster()) 00135 { 00136 p_file = output_file_handler.OpenOutputFile(stream.str()); 00137 } 00138 for (unsigned i = 0; i < mNumTimeSteps; i++) 00139 { 00140 double pseudo_ecg_at_one_timestep = ComputePseudoEcgAtOneTimeStep(i); 00141 if (PetscTools::AmMaster()) 00142 { 00143 *p_file << pseudo_ecg_at_one_timestep << "\n"; 00144 } 00145 } 00146 if (PetscTools::AmMaster()) 00147 { 00148 //write provenance info 00149 std::string comment = "# " + ChasteBuildInfo::GetProvenanceString(); 00150 *p_file << comment; 00151 p_file->close(); 00152 } 00153 } 00155 // Explicit instantiation 00157 00158 template class PseudoEcgCalculator<1,1,1>; 00159 template class PseudoEcgCalculator<1,2,1>; 00160 template class PseudoEcgCalculator<1,3,1>; 00161 template class PseudoEcgCalculator<2,2,1>; 00162 //template class PseudoEcgCalculator<2,3,1>; 00163 template class PseudoEcgCalculator<3,3,1>; 00164