ESyS-Particle  4.0.1
pi_storage_ed.hpp
00001 
00002 //                                                         //
00003 // Copyright (c) 2003-2011 by The University of Queensland //
00004 // Earth Systems Science Computational Centre (ESSCC)      //
00005 // http://www.uq.edu.au/esscc                              //
00006 //                                                         //
00007 // Primary Business: Brisbane, Queensland, Australia       //
00008 // Licensed under the Open Software License version 3.0    //
00009 // http://www.opensource.org/licenses/osl-3.0.php          //
00010 //                                                         //
00012 
00013 template<typename P,typename I>
00014 ParallelInteractionStorage_ED<P,I>::ParallelInteractionStorage_ED(AParallelParticleArray* ppa,const typename I::ParameterType& param):ParallelInteractionStorage_E<P,I>(ppa,param)
00015 {
00016   m_exIG=NULL;
00017   m_update_timestamp=0;
00018 }
00019 
00020 template<typename P,typename InteractionType>
00021 void  ParallelInteractionStorage_ED<P,InteractionType>::addExIG(AParallelInteractionStorage* eg)
00022 {
00023   console.Debug() << "setExIG " << eg << "\n";
00024   m_exIG=eg; 
00025   // clean out excluded interactions
00026   typename list<InteractionType>::iterator iter = this->m_interactions.begin();
00027   while(iter != this->m_interactions.end()){
00028     // check if in exIG
00029     vector<int> rm_pids=iter->getAllID();
00030     if(m_exIG->isIn(rm_pids)){
00031       console.XDebug() << "removing excluded: " << rm_pids[0] << " - " << rm_pids[1] << "\n";
00032       typename list<InteractionType>::iterator er_iter=iter;
00033       // get particle ids and remove pair from set
00034       this->m_set.erase(make_pair(rm_pids[0],rm_pids[1]));
00035       iter++;
00036       this->m_interactions.erase(er_iter);
00037     } else {
00038       iter++;
00039     }
00040   }
00041 }
00042 
00043 template<typename T,typename InteractionType>
00044 void ParallelInteractionStorage_ED<T,InteractionType>::setTimeStepSize(
00045   double dt
00046 )
00047 {
00048   console.Debug()
00049     << "setting time step size for "
00050     << this->m_interactions.size() << " interaction forces\n" ;
00051   this->m_param.setTimeStepSize(dt);
00052   for (
00053     typename std::list<InteractionType>::iterator it = this->m_interactions.begin();
00054     it != this->m_interactions.end();
00055     it++
00056   ){
00057     it->setTimeStepSize(dt);
00058   }
00059 }
00060 
00064 template<typename T,typename InteractionType>
00065 bool ParallelInteractionStorage_ED<T,InteractionType>::update()
00066 {
00067   //std::cout << "ParallelInteractionStorage_ED::update at node " << m_comm.rank() << std::endl << std::flush;
00068   int count_l=0;
00069   bool res=true;
00070 
00071   if (this->m_update_timestamp != this->m_ppa->getTimeStamp()){// m_ppa rebuild since last update 
00072     console.XDebug() << "node " << this->m_comm.rank() << " ppa has been rebuilt\n";
00073     // clean out old interactions if not flagged as persistent
00074     typename list<InteractionType>::iterator iter = this->m_interactions.begin();
00075     while(iter != this->m_interactions.end()){
00076       if(iter->isPersistent()){
00077         iter++;
00078         //console.XDebug() << "node " << m_comm.rank() << "persistent interaction\n";
00079       }else{
00080         typename list<InteractionType>::iterator er_iter=iter;
00081         // get particle ids and remove pair from set
00082         vector<int> rm_pids=iter->getAllID();
00083         this->m_set.erase(make_pair(rm_pids[0],rm_pids[1]));
00084         iter++;
00085         this->m_interactions.erase(er_iter);
00086       }
00087     }
00088     // get list  of pairs from m_ppa
00089     typename ParallelParticleArray<T>::PairListHandle plh =
00090       ((ParallelParticleArray<T>*)this->m_ppa)->getFullPairList();
00091     // generate interactions from pairs
00092     for(typename ParallelParticleArray<T>::PairListIterator iter=plh->begin();
00093         iter!=plh->end();
00094         iter++){
00095       // check vs. ExIG
00096       vector<int> tv;
00097       // ids in pair
00098       int id1=iter->first->getID();
00099       int id2=iter->second->getID();
00100       tv.push_back(id1);
00101       tv.push_back(id2);
00102       if(m_exIG!=NULL){ // if there is an ExIG
00103         if((!m_exIG->isIn(tv))&&(!this->isIn(tv))){  // if not already in or in ExIG
00104           this->m_interactions.push_back(
00105             InteractionType(iter->first,iter->second,this->m_param)
00106           );
00107           this->m_set.insert(make_pair(id1,id2));
00108           count_l++;
00109         }
00110       } else if (!(this->isIn(tv))) { // if no ExIG -> check only if alrady in
00111         this->m_interactions.push_back(
00112           InteractionType(iter->first,iter->second,this->m_param)
00113         );
00114         this->m_set.insert(make_pair(id1,id2));
00115       }
00116     }
00117   } else { // m_ppa not rebuild since last update -> just get additional interactions
00118     console.XDebug() << "node " << this->m_comm.rank() << " ppa not rebuilt\n";
00119     // get list  of pairs from m_ppa
00120     typename ParallelParticleArray<T>::PairListHandle plh =
00121       ((ParallelParticleArray<T>*)this->m_ppa)->getNewPairList();
00122     for (
00123       typename ParallelParticleArray<T>::PairListIterator iter=plh->begin();
00124       iter!=plh->end();
00125       iter++
00126     ){
00127       // check vs. ExIG
00128       vector<int> tv;
00129       // ids in pair
00130       int id1=iter->first->getID();
00131       int id2=iter->second->getID();
00132       tv.push_back(id1);
00133       tv.push_back(id2); 
00134       if(m_exIG!=NULL){
00135         if((!m_exIG->isIn(tv))&&(!(this->isIn(tv)))) {
00136           this->m_interactions.push_back(
00137             InteractionType(iter->first,iter->second, this->m_param)
00138           );
00139           this->m_set.insert(make_pair(id1,id2));
00140           count_l++;
00141         }
00142       } else if (!(this->isIn(tv))) {
00143         this->m_interactions.push_back(
00144           InteractionType(iter->first,iter->second,this->m_param)
00145         );
00146         this->m_set.insert(make_pair(id1,id2));
00147       }
00148     }
00149   }
00150   m_update_timestamp = this->m_ppa->getTimeStamp();
00151 
00152   console.Debug() << "added " << count_l << " pairs to ParallelInteractionStorage_ED\n";
00153   //  std::cout << "end ParallelInteractionStorage_ED::update at node " << m_comm.rank() << std::endl << std::flush;
00154 
00155         return res;
00156 }
00157 
00158 
00159 template<typename T,typename InteractionType>
00160 void ParallelInteractionStorage_ED<T,InteractionType>::calcHeatFrict()
00161 {
00162   console.Debug()
00163     << "calculating "
00164     << this->m_interactions.size()
00165     << " frictional heatings\n" ;
00166 
00167   for(
00168     typename list<InteractionType>::iterator it = this->m_interactions.begin();
00169     it != this->m_interactions.end();
00170     it++
00171   ){
00172     it->calcHeatFrict();
00173   }
00174 }
00175 
00176 template<typename P,typename InteractionType>
00177 void ParallelInteractionStorage_ED<P,InteractionType>::calcHeatTrans()
00178 {
00179   console.Debug()
00180     << "calculating "
00181     << this->m_interactions.size()
00182     << " heat transfers\n" ;
00183 
00184   for(
00185     typename list<InteractionType>::iterator it = this->m_interactions.begin();
00186     it != this->m_interactions.end();
00187     it++
00188   ){
00189     it->calcHeatTrans();
00190   }
00191 }
00192 
00196 template<typename P,typename InteractionType>
00197 void ParallelInteractionStorage_ED<P,InteractionType>::saveCheckPointData(std::ostream &oStream)
00198 {
00199   const std::string delim = "\n";
00200 
00201   typename ParallelInteractionStorage_E<P,InteractionType>::InteractionIterator it =
00202     this->getInnerInteractionIterator();
00203   oStream << InteractionType::getType() << delim;
00204   oStream << it.getNumRemaining();
00205   if (it.hasNext()) {
00206     oStream << delim;
00207     it.next().saveRestartData(oStream);
00208     while (it.hasNext())
00209     {
00210       oStream << delim;
00211       it.next().saveRestartData(oStream);
00212     }
00213   }
00214 }
00215 
00223 template<typename P,typename InteractionType>
00224 void ParallelInteractionStorage_ED<P,InteractionType>::loadCheckPointData(std::istream &iStream)
00225 {
00226   // read interaction type from stream
00227   std::string cp_interaction_type;
00228   iStream >> cp_interaction_type;
00229   // compare interaction type in stream with type of this IG
00230   // in not equal, signal error 
00231   if(cp_interaction_type!=InteractionType::getType()){
00232     std::cerr << "interaction types differ between checkpoint " 
00233               << cp_interaction_type << " and scipt "
00234               << InteractionType::getType() << std::endl;
00235   } else { // correct type -> read data
00236     // read nr. of bonds in IG
00237     int nconn;
00238     iStream >> nconn;
00239     std::cerr << "reading " << nconn << "  " << InteractionType::getType() << " interactions " << std::endl;
00240   
00241     // -- read bonds
00242     for(int i=0;i<nconn;i++){
00243       InteractionType new_bond;
00244       // read a bond
00245       new_bond.loadRestartData(iStream);
00246       // insert it into interaction storage
00247       this->tryInsert(new_bond);
00248     }
00249   }
00250 }