ESyS-Particle  4.0.1
trimesh_pis_eb.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 #include "Foundation/console.h"
00014 
00015 template<class ParticleType,class IType>
00016 const int TriMesh_PIS_EB<ParticleType,IType>::m_exchg_tag=44;
00017 
00025 template<class ParticleType,class IType>
00026 TriMesh_PIS_EB<ParticleType,IType>::TriMesh_PIS_EB(TriMesh* mesh_p,ParallelParticleArray<ParticleType>* ppa_p,typename IType::ParameterType param)
00027   : TriMesh_PIS<ParticleType>(mesh_p,ppa_p),m_comm(ppa_p->getComm())
00028 {
00029   console.XDebug() << "TriMesh_PIS_EB constructor\n";
00030   m_param=param;
00031   this->m_update_timestamp = 0;
00032 }
00033    
00034 template <class ParticleType,class IType> 
00035 void TriMesh_PIS_EB<ParticleType,IType>::setTimeStepSize(double dt)
00036 {
00037 }
00038 
00048 template <class ParticleType,class IType> 
00049 bool TriMesh_PIS_EB<ParticleType,IType>::isIn(const std::vector<int>& v)
00050 {
00051   bool res=false;
00052   
00053   if(v.size()<3){
00054     res=false;
00055   } else {
00056     switch (v[2]){
00057     case 0: res=m_tri_int_set.find(make_pair(v[0],v[1]))!=m_tri_int_set.end(); break;
00058     default: console.Error() << "wrong value in argument of TriMesh_PIS::isIn !!\n"; break;
00059    }
00060   }
00061 
00062   return res;
00063 }
00064 
00068 template<class ParticleType,class IType>
00069 void TriMesh_PIS_EB<ParticleType,IType>::calcForces()
00070 {
00071   console.XDebug() << "TriMesh_PIS_EB calculating " << m_triangle_interactions.size() << " triangle forces\n";
00072 
00073   // calculate forces for triangle interactions
00074   for(typename list<typename IType::TriIntType>::iterator tri_iter=m_triangle_interactions.begin();
00075       tri_iter!=m_triangle_interactions.end();
00076       tri_iter++){
00077     tri_iter->calcForces();
00078   }
00079 }
00080 
00083 template<class ParticleType,class IType>
00084 bool TriMesh_PIS_EB<ParticleType,IType>::update()
00085 {
00086   console.XDebug() << "TriMesh_PIS_EB::update on node " << m_comm.rank() << "\n"; 
00087   bool res=false;
00088 
00089   typename list<typename IType::TriIntType>::iterator iter=m_triangle_interactions.begin();
00090   while(iter!=m_triangle_interactions.end()){
00091     if(iter->broken()){
00092       res=true;
00093       typename list<typename IType::TriIntType>::iterator er_iter=iter;
00094       iter++;
00095       // remove ids from map
00096       m_tri_int_set.erase(make_pair(er_iter->getTid(),er_iter->getPid()));
00097       // remove interaction
00098       m_triangle_interactions.erase(er_iter);
00099     } else {
00100       iter++;
00101     }
00102   }
00103 
00104   console.XDebug() << "end TriMesh_PIS_EB::update on node " << m_comm.rank() << "\n"; 
00105   return res;
00106 }
00107    
00114 template<class ParticleType,class IType>
00115 void TriMesh_PIS_EB<ParticleType,IType>::exchange_boundary(int dim,int dir)
00116 {
00117   console.XDebug() << "TriMesh_PIS_EB::exchange_boundary(" << dim << "," << dir << ") at node " << m_comm.rank() << "\n";
00118   
00119   std::set<int> bdry_ids;
00120   std::vector<typename IType::TriIntType> recv_tri_buffer;
00121   std::vector<typename IType::TriIntType> send_tri_buffer;
00122 
00123   // --- setup data to send ---
00124   // get boundary
00125   bdry_ids = this->m_ppa->getBoundarySlabIds(dim,dir);
00126   // for all interactions
00127   for(typename list<typename IType::TriIntType>::iterator iter=m_triangle_interactions.begin();
00128       iter!=m_triangle_interactions.end();
00129       iter++){
00130     int pid=iter->getPid(); 
00131     if(bdry_ids.find(pid)!=bdry_ids.end()) { // if particle in interaction is in bdry -> put in buffer
00132       send_tri_buffer.push_back(*iter);
00133     }  
00134   }
00135   // ---- shift ----
00136   m_comm.shift_cont_packed(send_tri_buffer,recv_tri_buffer,dim,dir,m_exchg_tag);
00137   // insert received data
00138   for(typename std::vector<typename IType::TriIntType>::iterator iter=recv_tri_buffer.begin();
00139       iter!=recv_tri_buffer.end();
00140       iter++){
00141     tryInsert(*iter);
00142   }
00143   
00144   console.XDebug() << "end TriMesh_PIS_EB::exchange_boundary\n";
00145 }
00146 
00149 template<class ParticleType,class IType>
00150 void TriMesh_PIS_EB<ParticleType,IType>::exchange()
00151 {
00152   console.XDebug() << "TriMesh_PIS_EB::exchange\n";
00153   for(int i=0;i<3;i++){
00154     if(m_comm.get_dim(i)>1){
00155       // -- up --
00156       exchange_boundary(i,1);
00157       // -- down --
00158       exchange_boundary(i,-1);
00159     }
00160   }
00161   console.XDebug() << "end TriMesh_PIS_EB::exchange\n";
00162 }
00163    
00169 template<class ParticleType,class IType>
00170 void TriMesh_PIS_EB<ParticleType,IType>::rebuild()
00171 {
00172   console.XDebug() << "TriMesh_PIS_EB::rebuild on node " << m_comm.rank() << "\n"; 
00173   ParallelParticleArray<ParticleType>* t_ppa =
00174     (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be a dynamic_cast
00175   // for all triangle interactions
00176   typename list<typename IType::TriIntType>::iterator ti_iter=m_triangle_interactions.begin();
00177   while(ti_iter!=m_triangle_interactions.end()){
00178     int pid=ti_iter->getPid();
00179     ParticleType *part_p=t_ppa->getParticlePtrByIndex(pid);
00180     if(part_p!=NULL) { // particle is available in m_ppa -> setup pointer in interaction
00181       ti_iter->setPP(part_p);
00182       Triangle *tri_p = this->m_mesh->getTriangleById(ti_iter->getTid());
00183       ti_iter->setTP(tri_p);
00184       ti_iter++;
00185     } else { // particle is not available in m_ppa -> remove interaction
00186       const typename list<typename IType::TriIntType>::iterator er_iter=ti_iter;
00187       ti_iter++;
00188       m_tri_int_set.erase(make_pair(er_iter->getTid(),pid));
00189       m_triangle_interactions.erase(er_iter); 
00190     }
00191   }
00192 
00193   console.XDebug() << "end TriMesh_PIS_EB::rebuild on node " << m_comm.rank() << "\n"; 
00194 }
00195    
00207 template<class ParticleType,class IType>
00208 void TriMesh_PIS_EB<ParticleType,IType>::tryInsert(const std::vector<int>& pids)
00209 {
00210   ParallelParticleArray<ParticleType>* t_ppa =
00211     (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be dynamic_cast
00212 
00213   if(pids.size()<3){
00214     bool is_in=isIn(pids); // interaction already in
00215     ParticleType *part_p=t_ppa->getParticlePtrByIndex(pids[1]);
00216     if((!is_in) && (part_p!=NULL)){ 
00217       switch (pids[2]){
00218         case 0: {
00219           Triangle *tri_p = this->m_mesh->getTriangleById(pids[0]);
00220           if (tri_p!=NULL){
00221             m_triangle_interactions.push_back(
00222               typename IType::TriIntType(
00223                 part_p,
00224                 tri_p,
00225                 m_param,
00226                 this->m_ppa->isInInner(part_p->getPos())
00227               )
00228             );
00229             m_tri_int_set.insert(make_pair(pids[0],pids[1]));
00230           }
00231         } break;
00232         default : {
00233           console.Error()
00234             << "Error: pids[2]= " << pids[2]
00235             << "\n"; // Can't happen !! 
00236         }
00237       }
00238     }
00239   }
00240 }
00241 
00244 template<class ParticleType,class IType>
00245 void TriMesh_PIS_EB<ParticleType,IType>::tryInsert(const typename IType::TriIntType& In)
00246 {
00247   ParallelParticleArray<ParticleType>* t_ppa =
00248     (ParallelParticleArray<ParticleType>*)(this->m_ppa);
00249   // check if interaction is already in 
00250   bool is_in=(m_tri_int_set.find(make_pair(In.getTid(),In.getPid()))!=m_tri_int_set.end());
00251   ParticleType *part_p=t_ppa->getParticlePtrByIndex(In.getPid());
00252   if((!is_in) && (part_p!=NULL)){
00253     m_triangle_interactions.push_back(In);
00254     m_tri_int_set.insert(make_pair(In.getTid(),In.getPid()));
00255   }
00256 }
00257 
00260 template<class ParticleType,class IType>
00261 void TriMesh_PIS_EB<ParticleType,IType>::buildFromPPATagged(int tag,int mask)
00262 {
00263   console.XDebug() << "TriMesh_PIS_EB::buildFromPPATagged(" << tag << "," << mask << ")\n";
00264   
00265   console.XDebug() << "end  TriMesh_PIS_EB::buildFromPPATagged()";
00266 }
00267 
00270 template<class ParticleType,class IType>
00271 void TriMesh_PIS_EB<ParticleType,IType>::buildFromPPAByGap(double gmax)
00272 {
00273   console.XDebug() << "TriMesh_PIS_EB::buildFromPPAByGap(" << gmax << ")\n";
00274   set<int> id_set;
00275 
00276   // for all triangles
00277   for (
00278     TriMesh::triangle_iterator tri_iter = this->m_mesh->triangles_begin();
00279     tri_iter != this->m_mesh->triangles_end();
00280     tri_iter++
00281   ){
00282     // get particles near triangle
00283     typename ParallelParticleArray<ParticleType>::ParticleListHandle plh=
00284       ((ParallelParticleArray<ParticleType>*)this->m_ppa)->getParticlesNearTriangle(*tri_iter); 
00285     console.XDebug() << "triangle " << tri_iter->getID() << " nr. of particles : " << plh->size() << "\n"; 
00286     // for all particles found
00287     for(typename ParallelParticleArray<ParticleType>::ParticleListIterator p_iter=plh->begin();
00288         p_iter!=plh->end();
00289         p_iter++){
00290       // if valid interaction
00291       console.XDebug() << "interaction : " << tri_iter->getID() << " " << (*p_iter)->getID() << "\n";
00292       if(id_set.find((*p_iter)->getID())==id_set.end()){
00293         pair<bool,double> dist=tri_iter->dist((*p_iter)->getPos());
00294         console.XDebug() << "is valid: " << dist.first << " dist : " << dist.second << "\n";
00295         if(dist.first){
00296           // check gap
00297           double gap=fabs(dist.second-(*p_iter)->getRad());
00298           console.XDebug() << "radius: " << (*p_iter)->getRad() << " gap : " << gap << "\n";
00299           // if small enough, add
00300           if(gap<gmax){
00301             console.XDebug() << "Insert !!\n";
00302             bool in_flag = this->m_ppa->isInInner((*p_iter)->getPos());
00303             m_triangle_interactions.push_back(typename IType::TriIntType((*p_iter),&(*tri_iter),m_param,in_flag));
00304             m_tri_int_set.insert(make_pair(tri_iter->getID(),(*p_iter)->getID()));
00305             id_set.insert((*p_iter)->getID());
00306           }
00307         }
00308       }
00309     }
00310   }
00311   console.XDebug() << "end  TriMesh_PIS_EB::buildFromPPAByGap()";
00312 }
00313 
00319 template<class ParticleType,class IType>
00320 void TriMesh_PIS_EB<ParticleType,IType>::saveSnapShotData(std::ostream &oStream)
00321 {
00322   const std::string delim = "\n";
00323   typedef typename IType::TriIntType::CheckPointable CheckPointable;
00324 
00325   // stage 1: count how many interactions with "inner" particles we have
00326   int icount=0;
00327   for(typename list<typename IType::TriIntType>::iterator it=m_triangle_interactions.begin();
00328       it!=m_triangle_interactions.end();
00329       it++){
00330     if(it->isInner()) icount++;
00331   }
00332 
00333   // stage 2: write data
00334   oStream << IType::getType() << delim;
00335   oStream << icount << delim;
00336   for(typename list<typename IType::TriIntType>::iterator it=m_triangle_interactions.begin();
00337       it!=m_triangle_interactions.end();
00338       it++){
00339     if(it->isInner()) CheckPointable(*it).saveCheckPointData(oStream);
00340     oStream << delim;
00341   }
00342 }