PostProcessingWriter.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include "UblasCustomFunctions.hpp"
00030 #include "HeartConfig.hpp"
00031 #include "PostProcessingWriter.hpp"
00032 #include "PetscTools.hpp"
00033 #include "OutputFileHandler.hpp"
00034 #include "DistanceMapCalculator.hpp"
00035 #include "Version.hpp"
00036 
00037 #include <iostream>
00038 
00039 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00040 PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::PostProcessingWriter(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>& rMesh, std::string directory, std::string hdf5File, bool makeAbsolute)
00041                                 : mrMesh(rMesh)
00042 {
00043     mLo=mrMesh.GetDistributedVectorFactory()->GetLow();
00044     mHi=mrMesh.GetDistributedVectorFactory()->GetHigh();
00045     mpDataReader = new Hdf5DataReader(directory, hdf5File, makeAbsolute);
00046     mpCalculator = new PropagationPropertiesCalculator(mpDataReader);
00047     //check that the hdf file was generated by simulations from the same mesh
00048     assert(mpDataReader->GetNumberOfRows() == mrMesh.GetNumNodes());
00049 }
00050 
00051 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00052 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WritePostProcessingFiles()
00053 {
00054 // Please note that only the master processor should write to file.
00055 // Each of the private methods called here takes care of checking.
00056 
00057     if (HeartConfig::Instance()->IsApdMapsRequested())
00058     {
00059         std::vector<std::pair<double,double> > apd_maps;
00060         HeartConfig::Instance()->GetApdMaps(apd_maps);
00061         for (unsigned i=0; i<apd_maps.size(); i++)
00062         {
00063             WriteApdMapFile(apd_maps[i].first, apd_maps[i].second);
00064         }
00065     }
00066 
00067     if (HeartConfig::Instance()->IsUpstrokeTimeMapsRequested())
00068     {
00069         std::vector<double> upstroke_time_maps;
00070         HeartConfig::Instance()->GetUpstrokeTimeMaps(upstroke_time_maps);
00071         for (unsigned i=0; i<upstroke_time_maps.size(); i++)
00072         {
00073             WriteUpstrokeTimeMap(upstroke_time_maps[i]);
00074         }
00075     }
00076 
00077     if (HeartConfig::Instance()->IsMaxUpstrokeVelocityMapRequested())
00078     {
00079         std::vector<double> upstroke_velocity_maps;
00080         HeartConfig::Instance()->GetMaxUpstrokeVelocityMaps(upstroke_velocity_maps);
00081         for (unsigned i=0; i<upstroke_velocity_maps.size(); i++)
00082         {
00083             WriteMaxUpstrokeVelocityMap(upstroke_velocity_maps[i]);
00084         }
00085     }
00086 
00087     if (HeartConfig::Instance()->IsConductionVelocityMapsRequested())
00088     {
00089         std::vector<unsigned> conduction_velocity_maps;
00090         HeartConfig::Instance()->GetConductionVelocityMaps(conduction_velocity_maps);
00091 
00092         //get the mesh here
00093         DistanceMapCalculator<ELEMENT_DIM, SPACE_DIM> dist_map_calculator(mrMesh);
00094 
00095         for (unsigned i=0; i<conduction_velocity_maps.size(); i++)
00096         {
00097             std::vector<double> distance_map;
00098             std::vector<unsigned> origin_surface;
00099             origin_surface.push_back(conduction_velocity_maps[i]);
00100             dist_map_calculator.ComputeDistanceMap(origin_surface, distance_map);
00101             WriteConductionVelocityMap(conduction_velocity_maps[i], distance_map);
00102         }
00103     }
00104 }
00105 
00106 
00107 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00108 PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::~PostProcessingWriter()
00109 {
00110     delete mpDataReader;
00111     delete mpCalculator;
00112 }
00113 
00114 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00115 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteApdMapFile(double repolarisationPercentage, double threshold)
00116 {
00117     std::vector<std::vector<double> > output_data;
00118     //Fill in data
00119     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00120     {
00121         std::vector<double> apds;
00122         try
00123         {
00124             apds = mpCalculator->CalculateAllActionPotentialDurations(repolarisationPercentage, node_index, threshold);
00125             assert(apds.size() != 0);
00126         }
00127         catch(Exception& e)
00128         {
00129             assert(e.GetShortMessage()=="No full action potential was recorded" ||
00130                    e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage.");
00131             apds.push_back(0);
00132             assert(apds.size() == 1);
00133         }
00134         output_data.push_back(apds);
00135     }
00136     std::stringstream filename_stream;
00137     filename_stream << "Apd_" << repolarisationPercentage << "_" << threshold << "_Map.dat";
00138     WriteGenericFile(output_data, filename_stream.str());
00139 }
00140 
00141 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00142 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteUpstrokeTimeMap(double threshold)
00143 {
00144     std::vector<std::vector<double> > output_data;
00145     //Fill in data
00146     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00147     {
00148         std::vector<double> upstroke_times;
00149         try
00150         {
00151             upstroke_times = mpCalculator->CalculateUpstrokeTimes(node_index, threshold);
00152             assert(upstroke_times.size() != 0);
00153         }
00154         catch(Exception& e)
00155         {
00156             upstroke_times.push_back(0);
00157             assert(upstroke_times.size() == 1);
00158         }
00159         output_data.push_back(upstroke_times);
00160     }
00161     std::stringstream filename_stream;
00162     filename_stream << "UpstrokeTimeMap_" << threshold << ".dat";
00163     WriteGenericFile(output_data, filename_stream.str());
00164 }
00165 
00166 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00167 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteMaxUpstrokeVelocityMap(double threshold)
00168 {
00169     std::vector<std::vector<double> > output_data;
00170     //Fill in data
00171     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00172     {
00173         std::vector<double> upstroke_velocities;
00174         try
00175         {
00176             upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
00177             assert(upstroke_velocities.size() != 0);
00178         }
00179         catch(Exception& e)
00180         {
00181             upstroke_velocities.push_back(0);
00182             assert(upstroke_velocities.size() ==1);
00183         }
00184         output_data.push_back(upstroke_velocities);
00185     }
00186     std::stringstream filename_stream;
00187     filename_stream << "MaxUpstrokeVelocityMap_" << threshold << ".dat";
00188     WriteGenericFile(output_data, filename_stream.str());
00189 }
00190 
00191 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00192 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteConductionVelocityMap(unsigned originNode, std::vector<double> distancesFromOriginNode)
00193 {
00194     //Fill in data
00195     std::vector<std::vector<double> > output_data;
00196     for (unsigned dest_node = mLo; dest_node < mHi; dest_node++)
00197     {
00198         std::vector<double> conduction_velocities;
00199         try
00200         {
00201             conduction_velocities = mpCalculator->CalculateAllConductionVelocities(originNode, dest_node, distancesFromOriginNode[dest_node]);
00202             assert(conduction_velocities.size() != 0);
00203         }
00204         catch(Exception& e)
00205         {
00206             conduction_velocities.push_back(0);
00207             assert(conduction_velocities.size() == 1);
00208         }
00209         output_data.push_back(conduction_velocities);
00210     }
00211     std::stringstream filename_stream;
00212     filename_stream << "ConductionVelocityFromNode" << originNode << ".dat";
00213     WriteGenericFile(output_data, filename_stream.str());
00214 }
00215 
00216 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00217 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteAboveThresholdDepolarisationFile(double threshold )
00218 {
00219     std::vector<std::vector<double> > output_data;
00220 
00221     //Fill in data
00222     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00223     {
00224         std::vector<double> upstroke_velocities;
00225         std::vector<unsigned> above_threshold_depolarisations;
00226         std::vector<double> output_item;
00227         bool no_upstroke_occurred = false;
00228 
00229         try
00230         {
00231             upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
00232             assert(upstroke_velocities.size() != 0);
00233         }
00234         catch(Exception& e)
00235         {
00236             upstroke_velocities.push_back(0);
00237             assert(upstroke_velocities.size() ==1);
00238             no_upstroke_occurred = true;
00239         }
00240         //this method won't throw any exception, so there is no need to put it into the try/catch
00241         above_threshold_depolarisations =  mpCalculator->CalculateAllAboveThresholdDepolarisations(node_index, threshold);
00242 
00243         //count the total above threshold depolarisations
00244         unsigned total_number_of_above_threshold_depolarisations = 0;
00245         for (unsigned ead_index = 0; ead_index< above_threshold_depolarisations.size();ead_index++)
00246         {
00247             total_number_of_above_threshold_depolarisations = total_number_of_above_threshold_depolarisations + above_threshold_depolarisations[ead_index];
00248         }
00249 
00250         //for this item, puch back the number of upstrokes...
00251         if (no_upstroke_occurred)
00252         {
00253             output_item.push_back(0);
00254         }
00255         else
00256         {
00257             output_item.push_back(upstroke_velocities.size());
00258         }
00259         //... and the number of above thrshold depolarisations
00260         output_item.push_back((double) total_number_of_above_threshold_depolarisations);
00261 
00262         output_data.push_back(output_item);
00263     }
00264     std::stringstream filename_stream;
00265     filename_stream << "AboveThresholdDepolarisations" << threshold << ".dat";
00266     WriteGenericFile(output_data, filename_stream.str());
00267 }
00268 
00269 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00270 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteGenericFile(std::vector<std::vector<double> >& rDataPayload, std::string fileName)
00271 {
00272     OutputFileHandler output_file_handler(HeartConfig::Instance()->GetOutputDirectory() + "/output", false);
00273     for (unsigned writing_process=0; writing_process<PetscTools::GetNumProcs(); writing_process++)
00274     {
00275         if(PetscTools::GetMyRank() == writing_process)
00276         {
00277             out_stream p_file=out_stream(NULL);
00278             if (PetscTools::AmMaster())
00279             {
00280                 //Open the file for the first time
00281                 p_file = output_file_handler.OpenOutputFile(fileName);
00282                 //write provenance info
00283                 std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
00284                 *p_file << comment;
00285             }
00286             else
00287             {
00288                 //Append to the existing file
00289                 p_file = output_file_handler.OpenOutputFile(fileName, std::ios::app);
00290             }
00291             for (unsigned line_number=0; line_number<rDataPayload.size(); line_number++)
00292             {
00293                 for (unsigned i = 0; i < rDataPayload[line_number].size(); i++)
00294                 {
00295                     *p_file << rDataPayload[line_number][i] << "\t";
00296                 }
00297                 *p_file << std::endl;
00298             }
00299             p_file->close();
00300         }
00301         //Process i+1 waits for process i to close the file
00302         PetscTools::Barrier();
00303     }//Loop in writing_process
00304 }
00305 
00307 // Explicit instantiation
00309 
00310 template class PostProcessingWriter<1,1>;
00311 template class PostProcessingWriter<1,2>;
00312 template class PostProcessingWriter<2,2>;
00313 template class PostProcessingWriter<1,3>;
00314 //template class PostProcessingWriter<2,3>;
00315 template class PostProcessingWriter<3,3>;

Generated by  doxygen 1.6.2