ESyS-Particle
4.0.1
|
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 }