PostProcessingWriter.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
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 "PseudoEcgCalculator.hpp"
00036 #include "Version.hpp"
00037 
00038 #include <iostream>
00039 
00040 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00041 PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::PostProcessingWriter(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>& rMesh,
00042                                                                    std::string directory,
00043                                                                    std::string hdf5File,
00044                                                                    bool makeAbsolute,
00045                                                                    std::string voltageName)
00046         : mDirectory(directory),
00047           mHdf5File(hdf5File),
00048           mMakeAbsolute(makeAbsolute),
00049           mVoltageName(voltageName),
00050           mrMesh(rMesh)
00051 {
00052     mLo = mrMesh.GetDistributedVectorFactory()->GetLow();
00053     mHi = mrMesh.GetDistributedVectorFactory()->GetHigh();
00054     mpDataReader = new Hdf5DataReader(directory, hdf5File, makeAbsolute);
00055     mpCalculator = new PropagationPropertiesCalculator(mpDataReader,voltageName);
00056     //check that the hdf file was generated by simulations from the same mesh
00057     assert(mpDataReader->GetNumberOfRows() == mrMesh.GetNumNodes());
00058 }
00059 
00060 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00061 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WritePostProcessingFiles()
00062 {
00063 // Please note that only the master processor should write to file.
00064 // Each of the private methods called here takes care of checking.
00065 
00066     if (HeartConfig::Instance()->IsApdMapsRequested())
00067     {
00068         std::vector<std::pair<double,double> > apd_maps;
00069         HeartConfig::Instance()->GetApdMaps(apd_maps);
00070         for (unsigned i=0; i<apd_maps.size(); i++)
00071         {
00072             WriteApdMapFile(apd_maps[i].first, apd_maps[i].second);
00073         }
00074     }
00075 
00076     if (HeartConfig::Instance()->IsUpstrokeTimeMapsRequested())
00077     {
00078         std::vector<double> upstroke_time_maps;
00079         HeartConfig::Instance()->GetUpstrokeTimeMaps(upstroke_time_maps);
00080         for (unsigned i=0; i<upstroke_time_maps.size(); i++)
00081         {
00082             WriteUpstrokeTimeMap(upstroke_time_maps[i]);
00083         }
00084     }
00085 
00086     if (HeartConfig::Instance()->IsMaxUpstrokeVelocityMapRequested())
00087     {
00088         std::vector<double> upstroke_velocity_maps;
00089         HeartConfig::Instance()->GetMaxUpstrokeVelocityMaps(upstroke_velocity_maps);
00090         for (unsigned i=0; i<upstroke_velocity_maps.size(); i++)
00091         {
00092             WriteMaxUpstrokeVelocityMap(upstroke_velocity_maps[i]);
00093         }
00094     }
00095 
00096     if (HeartConfig::Instance()->IsConductionVelocityMapsRequested())
00097     {
00098         std::vector<unsigned> conduction_velocity_maps;
00099         HeartConfig::Instance()->GetConductionVelocityMaps(conduction_velocity_maps);
00100 
00101         //get the mesh here
00102         DistanceMapCalculator<ELEMENT_DIM, SPACE_DIM> dist_map_calculator(mrMesh);
00103 
00104         for (unsigned i=0; i<conduction_velocity_maps.size(); i++)
00105         {
00106             std::vector<double> distance_map;
00107             std::vector<unsigned> origin_surface;
00108             origin_surface.push_back(conduction_velocity_maps[i]);
00109             dist_map_calculator.ComputeDistanceMap(origin_surface, distance_map);
00110             WriteConductionVelocityMap(conduction_velocity_maps[i], distance_map);
00111         }
00112     }
00113     
00114     if (HeartConfig::Instance()->IsAnyNodalTimeTraceRequested())
00115     {
00116         std::vector<unsigned> requested_nodes;
00117         HeartConfig::Instance()->GetNodalTimeTraceRequested(requested_nodes);
00118         WriteVariablesOverTimeAtNodes(requested_nodes);
00119     }
00120 
00121     if (HeartConfig::Instance()->IsPseudoEcgCalculationRequested())
00122     {
00123         std::vector<ChastePoint<SPACE_DIM> > electrodes;
00124         HeartConfig::Instance()->GetPseudoEcgElectrodePositions(electrodes);
00125         
00126         for (unsigned i=0; i<electrodes.size(); i++)
00127         {
00128             PseudoEcgCalculator<ELEMENT_DIM,SPACE_DIM,1> calculator(mrMesh, electrodes[i], mDirectory, mHdf5File, mVoltageName, mMakeAbsolute);
00129             calculator.WritePseudoEcg();
00130         }
00131     }
00132 }
00133 
00134 
00135 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00136 PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::~PostProcessingWriter()
00137 {
00138     delete mpDataReader;
00139     delete mpCalculator;
00140 }
00141 
00142 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00143 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteApdMapFile(double repolarisationPercentage, double threshold)
00144 {
00145     std::vector<std::vector<double> > output_data;
00146     //Fill in data
00147     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00148     {
00149         std::vector<double> apds;
00150         try
00151         {
00152             apds = mpCalculator->CalculateAllActionPotentialDurations(repolarisationPercentage, node_index, threshold);
00153             assert(apds.size() != 0);
00154         }
00155         catch(Exception& e)
00156         {
00157             assert(e.GetShortMessage()=="No full action potential was recorded" ||
00158                    e.GetShortMessage()=="AP did not occur, never exceeded threshold voltage.");
00159             apds.push_back(0);
00160             assert(apds.size() == 1);
00161         }
00162         output_data.push_back(apds);
00163     }
00164     std::stringstream filename_stream;
00165     filename_stream << "Apd_" << repolarisationPercentage << "_" << threshold << "_Map.dat";
00166     WriteGenericFile(output_data, filename_stream.str());
00167 }
00168 
00169 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00170 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteUpstrokeTimeMap(double threshold)
00171 {
00172     std::vector<std::vector<double> > output_data;
00173     //Fill in data
00174     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00175     {
00176         std::vector<double> upstroke_times;
00177         try
00178         {
00179             upstroke_times = mpCalculator->CalculateUpstrokeTimes(node_index, threshold);
00180             assert(upstroke_times.size() != 0);
00181         }
00182         catch(Exception& e)
00183         {
00184             upstroke_times.push_back(0);
00185             assert(upstroke_times.size() == 1);
00186         }
00187         output_data.push_back(upstroke_times);
00188     }
00189     std::stringstream filename_stream;
00190     filename_stream << "UpstrokeTimeMap_" << threshold << ".dat";
00191     WriteGenericFile(output_data, filename_stream.str());
00192 }
00193 
00194 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00195 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteMaxUpstrokeVelocityMap(double threshold)
00196 {
00197     std::vector<std::vector<double> > output_data;
00198     //Fill in data
00199     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00200     {
00201         std::vector<double> upstroke_velocities;
00202         try
00203         {
00204             upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
00205             assert(upstroke_velocities.size() != 0);
00206         }
00207         catch(Exception& e)
00208         {
00209             upstroke_velocities.push_back(0);
00210             assert(upstroke_velocities.size() ==1);
00211         }
00212         output_data.push_back(upstroke_velocities);
00213     }
00214     std::stringstream filename_stream;
00215     filename_stream << "MaxUpstrokeVelocityMap_" << threshold << ".dat";
00216     WriteGenericFile(output_data, filename_stream.str());
00217 }
00218 
00219 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00220 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteConductionVelocityMap(unsigned originNode, std::vector<double> distancesFromOriginNode)
00221 {
00222     //Fill in data
00223     std::vector<std::vector<double> > output_data;
00224     for (unsigned dest_node = mLo; dest_node < mHi; dest_node++)
00225     {
00226         std::vector<double> conduction_velocities;
00227         try
00228         {
00229             conduction_velocities = mpCalculator->CalculateAllConductionVelocities(originNode, dest_node, distancesFromOriginNode[dest_node]);
00230             assert(conduction_velocities.size() != 0);
00231         }
00232         catch(Exception& e)
00233         {
00234             conduction_velocities.push_back(0);
00235             assert(conduction_velocities.size() == 1);
00236         }
00237         output_data.push_back(conduction_velocities);
00238     }
00239     std::stringstream filename_stream;
00240     filename_stream << "ConductionVelocityFromNode" << originNode << ".dat";
00241     WriteGenericFile(output_data, filename_stream.str());
00242 }
00243 
00244 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00245 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteAboveThresholdDepolarisationFile(double threshold )
00246 {
00247     std::vector<std::vector<double> > output_data;
00248 
00249     //Fill in data
00250     for (unsigned node_index = mLo; node_index < mHi; node_index++)
00251     {
00252         std::vector<double> upstroke_velocities;
00253         std::vector<unsigned> above_threshold_depolarisations;
00254         std::vector<double> output_item;
00255         bool no_upstroke_occurred = false;
00256 
00257         try
00258         {
00259             upstroke_velocities = mpCalculator->CalculateAllMaximumUpstrokeVelocities(node_index, threshold);
00260             assert(upstroke_velocities.size() != 0);
00261         }
00262         catch(Exception& e)
00263         {
00264             upstroke_velocities.push_back(0);
00265             assert(upstroke_velocities.size() ==1);
00266             no_upstroke_occurred = true;
00267         }
00268         //this method won't throw any exception, so there is no need to put it into the try/catch
00269         above_threshold_depolarisations =  mpCalculator->CalculateAllAboveThresholdDepolarisations(node_index, threshold);
00270 
00271         //count the total above threshold depolarisations
00272         unsigned total_number_of_above_threshold_depolarisations = 0;
00273         for (unsigned ead_index = 0; ead_index< above_threshold_depolarisations.size();ead_index++)
00274         {
00275             total_number_of_above_threshold_depolarisations = total_number_of_above_threshold_depolarisations + above_threshold_depolarisations[ead_index];
00276         }
00277 
00278         //for this item, puch back the number of upstrokes...
00279         if (no_upstroke_occurred)
00280         {
00281             output_item.push_back(0);
00282         }
00283         else
00284         {
00285             output_item.push_back(upstroke_velocities.size());
00286         }
00287         //... and the number of above thrshold depolarisations
00288         output_item.push_back((double) total_number_of_above_threshold_depolarisations);
00289 
00290         output_data.push_back(output_item);
00291     }
00292     std::stringstream filename_stream;
00293     filename_stream << "AboveThresholdDepolarisations" << threshold << ".dat";
00294     WriteGenericFile(output_data, filename_stream.str());
00295 }
00296 
00297 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00298 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteVariablesOverTimeAtNodes(std::vector<unsigned>& rNodeIndices)
00299 {
00300     std::vector<std::string> variable_names = mpDataReader->GetVariableNames();
00301 
00302     //we will write one file per variable in the hdf5 file
00303     for (unsigned name_index=0; name_index < variable_names.size(); name_index++)
00304     {
00305         std::vector<std::vector<double> > output_data;
00306         if (PetscTools::AmMaster())//only master process fills the data structure
00307         {
00308             //allocate memory: NXM matrix where N = numbe rof time stpes and M number of requested nodes
00309             output_data.resize( mpDataReader->GetUnlimitedDimensionValues().size() );
00310             for (unsigned j = 0; j < mpDataReader->GetUnlimitedDimensionValues().size(); j++)
00311             {
00312                 output_data[j].resize(rNodeIndices.size());
00313             }
00314 
00315             for (unsigned requested_index = 0; requested_index < rNodeIndices.size(); requested_index++)
00316             {
00317                 unsigned node_index = rNodeIndices[requested_index];
00318 
00319                 //handle permutation, if any
00320                 if ( (mrMesh.rGetNodePermutation().size() != 0) &&
00321                       !HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering() )
00322                 {
00323                     node_index = mrMesh.rGetNodePermutation()[ rNodeIndices[requested_index] ];
00324                 }
00325 
00326                 //grab the data from the hdf5 file.
00327                 std::vector<double> time_series =  mpDataReader->GetVariableOverTime(variable_names[name_index], node_index);
00328                 assert ( time_series.size() == mpDataReader->GetUnlimitedDimensionValues().size());
00329 
00330                 //fill the output_data data structure
00331                 for (unsigned time_step = 0; time_step < time_series.size(); time_step++)
00332                 {
00333                     output_data[time_step][requested_index] = time_series[time_step];
00334                 }
00335             }
00336         }
00337         std::stringstream filename_stream;
00338         filename_stream << "NodalTraces_" << variable_names[name_index] << ".dat";
00339         WriteGenericFile(output_data, filename_stream.str());
00340     }
00341 }
00342 
00343 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00344 void PostProcessingWriter<ELEMENT_DIM, SPACE_DIM>::WriteGenericFile(std::vector<std::vector<double> >& rDataPayload, std::string fileName)
00345 {
00346     OutputFileHandler output_file_handler(HeartConfig::Instance()->GetOutputDirectory() + "/output", false);
00347     for (unsigned writing_process=0; writing_process<PetscTools::GetNumProcs(); writing_process++)
00348     {
00349         if(PetscTools::GetMyRank() == writing_process)
00350         {
00351             out_stream p_file=out_stream(NULL);
00352             if (PetscTools::AmMaster())
00353             {
00354                 //Open the file for the first time
00355                 p_file = output_file_handler.OpenOutputFile(fileName);
00356                 //write provenance info
00357                 std::string comment = "# " + ChasteBuildInfo::GetProvenanceString();
00358                 *p_file << comment;
00359             }
00360             else
00361             {
00362                 //Append to the existing file
00363                 p_file = output_file_handler.OpenOutputFile(fileName, std::ios::app);
00364             }
00365             for (unsigned line_number=0; line_number<rDataPayload.size(); line_number++)
00366             {
00367                 for (unsigned i = 0; i < rDataPayload[line_number].size(); i++)
00368                 {
00369                     *p_file << rDataPayload[line_number][i] << "\t";
00370                 }
00371                 *p_file << std::endl;
00372             }
00373             p_file->close();
00374         }
00375         //Process i+1 waits for process i to close the file
00376         PetscTools::Barrier();
00377     }//Loop in writing_process
00378 }
00379 
00381 // Explicit instantiation
00383 
00384 template class PostProcessingWriter<1,1>;
00385 template class PostProcessingWriter<1,2>;
00386 template class PostProcessingWriter<2,2>;
00387 template class PostProcessingWriter<1,3>;
00388 //template class PostProcessingWriter<2,3>;
00389 template class PostProcessingWriter<3,3>;

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5