ESyS-Particle  4.0.1
mesh2d_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 template<class ParticleType,class IType>
00014 const int Mesh2D_PIS_EB<ParticleType,IType>::m_edge_exchg_tag=45;
00015 template<class ParticleType,class IType>
00016 const int Mesh2D_PIS_EB<ParticleType,IType>::m_corner_exchg_tag=46;
00017 
00025 template<class ParticleType,class IType>
00026 Mesh2D_PIS_EB<ParticleType,IType>::Mesh2D_PIS_EB(Mesh2D* mesh_p,ParallelParticleArray<ParticleType>* ppa_p,typename IType::ParameterType param)
00027   :Mesh2D_PIS<ParticleType>(mesh_p,ppa_p),m_comm(ppa_p->getComm())
00028 {
00029   console.XDebug() << "Mesh2D_PIS_EB constructor\n";
00030   m_param=param;
00031   this->m_update_timestamp=0;
00032 }
00033    
00043 template<class ParticleType,class IType>
00044 bool Mesh2D_PIS_EB<ParticleType,IType>::isIn(const std::vector<int>& v)
00045 {
00046   bool res=false;
00047   
00048   if(v.size()<3){
00049     res=false;
00050   } else {
00051     switch (v[2]){
00052     case 1: res=m_edge_int_set.find(make_pair(v[0],v[1]))!=m_edge_int_set.end(); break;
00053     case 2: res=m_edge_int_set.find(make_pair(v[0],v[1]))!=m_corner_int_set.end(); break;
00054     default: console.Error() << "wrong value in argument of Mesh2D_PIS::isIn !!\n"; break;
00055     }
00056   }
00057 
00058   return res;
00059 }
00060 
00064 template<class ParticleType,class IType>
00065 void Mesh2D_PIS_EB<ParticleType,IType>::calcForces()
00066 {
00067   console.XDebug() << "Mesh2D_PIS_EB calculating " << m_edge_interactions.size() << " edge forces and"
00068                    << m_corner_interactions.size() << " corner forces\n";
00069 
00070   // calculate forces for edge interactions
00071   for(typename list<typename IType::TriIntType>::iterator ed_iter=m_edge_interactions.begin();
00072       ed_iter!=m_edge_interactions.end();
00073       ed_iter++){
00074     ed_iter->calcForces();
00075   }
00076   // calculate forces for corner interactions
00077   for(typename list<typename IType::CornerIntType>::iterator c_iter=m_corner_interactions.begin();
00078       c_iter!=m_corner_interactions.end();
00079       c_iter++){
00080     c_iter->calcForces();
00081   }
00082 }
00083 
00086 template<class ParticleType,class IType>
00087 bool Mesh2D_PIS_EB<ParticleType,IType>::update()
00088 {
00089   console.XDebug() << "Mesh2D_PIS_EB::update on node " << m_comm.rank() << "\n"; 
00090   bool res=false;
00091 
00092   // edge interactions
00093   typename list<typename IType::TriIntType>::iterator iter=m_edge_interactions.begin();
00094   while(iter!=m_edge_interactions.end()){
00095     if(iter->broken()){
00096       res=true;
00097       typename list<typename IType::TriIntType>::iterator er_iter=iter;
00098       iter++;
00099       // remove ids from map
00100       m_edge_int_set.erase(make_pair(er_iter->getTid(),er_iter->getPid()));
00101       // remove interaction
00102       m_edge_interactions.erase(er_iter);
00103     } else {
00104       iter++;
00105     }
00106   }
00107   // corner interactions
00108   typename list<typename IType::CornerIntType>::iterator c_iter=m_corner_interactions.begin();
00109   while(c_iter!=m_corner_interactions.end()){
00110     if(c_iter->broken()){
00111       res=true;
00112       typename list<typename IType::CornerIntType>::iterator cr_iter=c_iter;
00113       c_iter++;
00114       // remove ids from map
00115       m_corner_int_set.erase(make_pair(cr_iter->getCid(),cr_iter->getPid()));
00116       // remove interaction
00117       m_corner_interactions.erase(cr_iter);
00118     } else {
00119       c_iter++;
00120     }
00121   }
00122   console.XDebug() << "end Mesh2D_PIS_EB::update on node " << m_comm.rank() << "\n"; 
00123   return res;
00124 }
00125 
00126 
00129 template<class ParticleType,class IType>
00130 void Mesh2D_PIS_EB<ParticleType,IType>::exchange()
00131 {
00132   console.XDebug() << "TriMesh_PIS_EB::exchange\n";
00133   for(int i=0;i<3;i++){
00134     if(m_comm.get_dim(i)>1){
00135       // -- up --
00136       exchange_boundary(i,1);
00137       // -- down --
00138       exchange_boundary(i,-1);
00139     }
00140   }
00141   console.XDebug() << "end TriMesh_PIS_EB::exchange\n";
00142 }
00143    
00150 template<class ParticleType,class IType>
00151 void Mesh2D_PIS_EB<ParticleType,IType>::exchange_boundary(int dim,int dir)
00152 {
00153   console.XDebug() << "Mesh2D_PIS_EB::exchange_boundary(" << dim << "," << dir << ") at node " << m_comm.rank() << "\n";
00154   
00155   std::set<int> bdry_ids;
00156   std::vector<typename IType::TriIntType> recv_tri_buffer;
00157   std::vector<typename IType::TriIntType> send_tri_buffer;
00158   std::vector<typename IType::CornerIntType> recv_corner_buffer;
00159   std::vector<typename IType::CornerIntType> send_corner_buffer;
00160 
00161   // --- setup data to send ---
00162   // get boundary
00163   bdry_ids = this->m_ppa->getBoundarySlabIds(dim,dir);
00164   // --- edges ---
00165   // for all interactions
00166   for(typename std::list<typename IType::TriIntType>::iterator iter=m_edge_interactions.begin();
00167       iter!=m_edge_interactions.end();
00168       iter++){
00169     int pid=iter->getPid(); 
00170     if(bdry_ids.find(pid)!=bdry_ids.end()) { // if particle in interaction is in bdry -> put in buffer
00171       send_tri_buffer.push_back(*iter);
00172     }  
00173   }
00174   // shift 
00175   m_comm.shift_cont_packed(send_tri_buffer,recv_tri_buffer,dim,dir,m_edge_exchg_tag);
00176   // insert received data
00177   for(typename std::vector<typename IType::TriIntType>::iterator iter=recv_tri_buffer.begin();
00178       iter!=recv_tri_buffer.end();
00179       iter++){
00180     tryInsert(*iter);
00181   }
00182   // --- corners ---
00183   // for all interactions
00184   for(typename std::list<typename IType::CornerIntType>::iterator iter=m_corner_interactions.begin();
00185        iter!=m_corner_interactions.end();
00186        iter++){
00187     int pid=iter->getPid(); 
00188     if(bdry_ids.find(pid)!=bdry_ids.end()) { // if particle in interaction is in bdry -> put in buffer
00189       send_corner_buffer.push_back(*iter);
00190     }  
00191   }
00192   // shift 
00193   m_comm.shift_cont_packed(send_corner_buffer,recv_corner_buffer,dim,dir,m_corner_exchg_tag);
00194   // insert received data
00195   for(typename std::vector<typename IType::CornerIntType>::iterator iter=recv_corner_buffer.begin();
00196       iter!=recv_corner_buffer.end();
00197       iter++){
00198     tryInsert(*iter);
00199   }
00200   
00201   console.XDebug() << "end Mesh2D_PIS_EB::exchange_boundary\n";
00202 }
00203 
00204 template<class ParticleType,class IType>
00205 void Mesh2D_PIS_EB<ParticleType,IType>::setTimeStepSize(double dt)
00206 {
00207 }
00208 
00214 template<class ParticleType,class IType>
00215 void Mesh2D_PIS_EB<ParticleType,IType>::rebuild()
00216 {
00217   console.XDebug() << "Mesh2D_PIS_EB::rebuild on node " << m_comm.rank() << "\n"; 
00218   ParallelParticleArray<ParticleType>* t_ppa =
00219     (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be a dynamic_cast
00220   // for all edge interactions
00221   typename std::list<typename IType::TriIntType>::iterator ti_iter=m_edge_interactions.begin();
00222   while(ti_iter!=m_edge_interactions.end()){
00223     int pid=ti_iter->getPid();
00224     ParticleType *part_p=t_ppa->getParticlePtrByIndex(pid);
00225     if(part_p!=NULL) { // particle is available in m_ppa -> setup pointer in interaction
00226       ti_iter->setPP(part_p);
00227       Edge2D *ed_p = this->m_mesh->getEdgeById(ti_iter->getTid());
00228       ti_iter->setTP(ed_p);
00229       ti_iter++;
00230     } else { // particle is not available in m_ppa -> remove interaction
00231       const typename list<typename IType::TriIntType>::iterator er_iter=ti_iter;
00232       ti_iter++;
00233       m_edge_int_set.erase(make_pair(er_iter->getTid(),pid));
00234       m_edge_interactions.erase(er_iter); 
00235     }
00236   }
00237   // and now for the corners
00238   typename list<typename IType::CornerIntType>::iterator ci_iter=m_corner_interactions.begin();
00239   while(ci_iter!=m_corner_interactions.end()){
00240     int pid=ci_iter->getPid();
00241     ParticleType *part_p=t_ppa->getParticlePtrByIndex(pid);
00242     if(part_p!=NULL) { // particle is available in m_ppa -> setup pointer in interaction
00243       ci_iter->setPP(part_p);
00244       Corner2D *co_p = this->m_mesh->getCornerById(ci_iter->getCid());
00245       ci_iter->setCP(co_p);
00246       ci_iter++;
00247     } else { // particle is not available in m_ppa -> remove interaction
00248       const typename list<typename IType::CornerIntType>::iterator cr_iter=ci_iter;
00249       ci_iter++;
00250       m_corner_int_set.erase(make_pair(cr_iter->getCid(),pid));
00251       m_corner_interactions.erase(cr_iter); 
00252     }
00253   }
00254 
00255   console.XDebug() << "end Mesh2D_PIS_EB::rebuild on node " << m_comm.rank() << "\n"; 
00256 }
00257    
00260 template<class ParticleType,class IType>
00261 void Mesh2D_PIS_EB<ParticleType,IType>::tryInsert(const typename IType::TriIntType& In)
00262 {
00263   console.XDebug() << "Mesh2D_PIS_EB::tryInsert(const typename IType::TriIntType& In)\n";
00264   ParallelParticleArray<ParticleType>* t_ppa =
00265     (ParallelParticleArray<ParticleType>*)(this->m_ppa);
00266   // check if interaction is already in 
00267   bool is_in=(m_edge_int_set.find(make_pair(In.getTid(),In.getPid()))!=m_edge_int_set.end());
00268   ParticleType *part_p=t_ppa->getParticlePtrByIndex(In.getPid());
00269   if((!is_in) && (part_p!=NULL)){
00270     m_edge_interactions.push_back(In);
00271     m_edge_int_set.insert(make_pair(In.getTid(),In.getPid()));
00272   }
00273 }
00274 
00277 template<class ParticleType,class IType>
00278 void Mesh2D_PIS_EB<ParticleType,IType>::tryInsert(const typename IType::CornerIntType& In)
00279 {
00280   console.XDebug() << "Mesh2D_PIS_EB::tryInsert(const typename IType::CornerIntType& In)\n";
00281   ParallelParticleArray<ParticleType>* t_ppa =
00282     (ParallelParticleArray<ParticleType>*)(this->m_ppa);
00283   // check if interaction is already in 
00284   bool is_in=(m_corner_int_set.find(make_pair(In.getCid(),In.getPid()))!=m_corner_int_set.end());
00285   ParticleType *part_p=t_ppa->getParticlePtrByIndex(In.getPid());
00286   if((!is_in) && (part_p!=NULL)){
00287     m_corner_interactions.push_back(In);
00288     m_corner_int_set.insert(make_pair(In.getCid(),In.getPid()));
00289   }
00290 }
00291 
00303 template<class ParticleType,class IType>
00304 void Mesh2D_PIS_EB<ParticleType,IType>::tryInsert(const vector<int>& pids)
00305 {
00306   console.XDebug() << "Mesh2D_PIS_EB::(const vector<int>& pids)\n";
00307   ParallelParticleArray<ParticleType>* t_ppa =
00308     (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be dynamic_cast
00309 
00310   if(pids.size()<3){
00311     bool is_in=isIn(pids); // interaction already in
00312     ParticleType *part_p=t_ppa->getParticlePtrByIndex(pids[1]);
00313     if((!is_in) && (part_p!=NULL)){ 
00314       switch (pids[2]){
00315       case 1: { // edge
00316         Edge2D *edge_p = this->m_mesh->getEdgeById(pids[0]);
00317         if(edge_p!=NULL){
00318           m_edge_interactions.push_back(
00319             typename IType::TriIntType(
00320               part_p,
00321               edge_p,
00322               m_param,
00323               this->m_ppa->isInInner(part_p->getPos())
00324             )
00325           );
00326           m_edge_int_set.insert(make_pair(pids[0],pids[1]));
00327         } else {
00328           console.Error() << "ERROR: Wrong edge id " << pids[0] << " in Mesh2D_PIS_EB::tryInsert\n";
00329         }
00330       } break;
00331       case 2: {
00332         Corner2D *corner_p = this->m_mesh->getCornerById(pids[0]);
00333         if(corner_p!=NULL){
00334           m_corner_interactions.push_back(
00335             typename IType::CornerIntType(
00336               part_p,
00337               corner_p,
00338               m_param,
00339               this->m_ppa->isInInner(part_p->getPos())
00340             )
00341           );
00342           m_corner_int_set.insert(make_pair(pids[0],pids[1]));
00343         } else {
00344           console.Error() << "ERROR: Wrong corner id " << pids[0] << " in Mesh2D_PIS_EB::tryInsert\n";
00345         }
00346       } break;
00347       default : console.Error() << "Error: pids[2]= " << pids[2] << "\n"; // Can't happen !! 
00348       }
00349     }
00350   }
00351 }
00352 
00353 
00356 template<class ParticleType,class IType>
00357 void Mesh2D_PIS_EB<ParticleType,IType>::buildFromPPATagged(int tag ,int mask)
00358 {
00359   console.XDebug() << "Mesh2D_PIS_EB::buildFromPPATagged(" << tag << "," << mask << ") not implemented \n";
00360   
00361   console.XDebug() << "end  Mesh2D_PIS_EB::buildFromPPATagged()";
00362 }
00363 
00369 template<class ParticleType,class IType>
00370 void Mesh2D_PIS_EB<ParticleType,IType>::buildFromPPAByGap(double gmax)
00371 {
00372   console.XDebug() << "Mesh2D_PIS_EB::buildFromPPAByGap(" << gmax << ")\n";
00373   set<int> id_set;
00374 
00375   // for all edges
00376   for(
00377     Mesh2D::edge_iterator ed_iter = this->m_mesh->edges_begin();
00378     ed_iter != this->m_mesh->edges_end();
00379     ed_iter++
00380   ){
00381     // get particles near edge
00382     typename ParallelParticleArray<ParticleType>::ParticleListHandle plh=
00383       ((ParallelParticleArray<ParticleType>*)this->m_ppa)->getParticlesNearEdge(&(*ed_iter)); 
00384     console.Debug() << "edge " << ed_iter->getID() << " nr. of particles : " << plh->size() << "\n"; 
00385     // for all particles found
00386     for(
00387       typename ParallelParticleArray<ParticleType>::ParticleListIterator p_iter=plh->begin();
00388       p_iter!=plh->end();
00389       p_iter++){
00390       // if valid interaction
00391       console.Debug() << "interaction : " << ed_iter->getID() << " " << (*p_iter)->getID() << "\n";
00392       if(id_set.find((*p_iter)->getID())==id_set.end()){
00393         pair<bool,double> dist=ed_iter->dist((*p_iter)->getPos());
00394         console.Debug() << "is valid: " << dist.first << " dist : " << dist.second << "\n";
00395         if(dist.first){
00396           // check ga
00397           double gap=fabs(dist.second-(*p_iter)->getRad());
00398           console.Debug() << "radius: " << (*p_iter)->getRad() << " gap : " << gap << "\n";
00399           // if small enough, add
00400           if(gap<gmax){
00401             console.Debug() << "Inserting " << (*p_iter)->getID() << " !!\n";
00402             bool in_flag = this->m_ppa->isInInner((*p_iter)->getPos());
00403             m_edge_interactions.push_back(typename IType::TriIntType((*p_iter),&(*ed_iter),m_param,in_flag));
00404             m_edge_int_set.insert(make_pair(ed_iter->getID(),(*p_iter)->getID()));
00405             id_set.insert((*p_iter)->getID());
00406           }
00407         }
00408       }
00409     }
00410   }
00411 }
00412 
00413 template <class ParticleType,class IType> 
00414 typename Mesh2D_PIS_EB<ParticleType,IType>::InteractionIterator  Mesh2D_PIS_EB<ParticleType,IType>::getInnerInteractionIterator()
00415 {
00416   return 
00417     InteractionIterator(
00418       m_edge_interactions.begin(),
00419       m_edge_interactions.end(),
00420       this->m_ppa
00421   );
00422 }
00423 
00424 template <class ParticleType,class IType> 
00425 void Mesh2D_PIS_EB<ParticleType,IType>::saveCheckPointData(std::ostream& ost)
00426 {
00427   const std::string delim = "\n";
00428   typedef typename IType::TriIntType::CheckPointable CheckPointable;
00429 
00430   InteractionIterator it = getInnerInteractionIterator();
00431   ost << IType::getType() << delim;
00432   ost << it.getNumRemaining();
00433   while (it.hasNext()){
00434     ost << delim;
00435     CheckPointable(it.next()).saveCheckPointData(ost);
00436   }
00437 }
00438 
00439 template <class ParticleType,class IType> 
00440 void Mesh2D_PIS_EB<ParticleType,IType>::loadCheckPointData(std::istream& ost)
00441 {
00442   console.Critical() << "Mesh2D_PIS_EB<ParticleType,IType>::loadCheckPointData NOT IMPLEMENTED\n";
00443 }
00444 
00445 // --- Iterator members ---
00446 
00447 template <class ParticleType,class IType> 
00448 Mesh2D_PIS_EB<ParticleType,IType>::InteractionIterator::InteractionIterator(Iterator begin, Iterator end, AParallelParticleArray* ppa)
00449   : m_numRemaining(0),
00450     m_it(end),
00451     m_end(end),
00452     m_ppa(ppa)
00453 {
00454   m_numRemaining = 0;
00455   for (Iterator it = begin; it != end; it++) {
00456     if  (isInner(it)) {
00457       m_numRemaining++;
00458     }
00459   }
00460   m_it  = begin;
00461   m_end = end;
00462 }
00463 
00464 
00465 template <class ParticleType,class IType> 
00466 bool Mesh2D_PIS_EB<ParticleType,IType>::InteractionIterator::hasNext()
00467 {
00468   return (m_numRemaining > 0);
00469 }
00470 
00471 template <class ParticleType,class IType> 
00472 typename IType::TriIntType& Mesh2D_PIS_EB<ParticleType,IType>::InteractionIterator::next()
00473 {
00474   while (!isInner(m_it)) {
00475     m_it++;
00476   }
00477   Interaction &i = *m_it;
00478   m_it++;
00479   m_numRemaining--;
00480   return i;
00481 }
00482 
00483 template <class ParticleType,class IType> 
00484 int Mesh2D_PIS_EB<ParticleType,IType>::InteractionIterator::getNumRemaining()
00485 {
00486   return m_numRemaining;
00487 }
00488 
00489 template <class ParticleType,class IType> 
00490 bool Mesh2D_PIS_EB<ParticleType,IType>::InteractionIterator::isInner(const Iterator& it)
00491 {
00492   return m_ppa->isInInner(it->getPos());
00493 }