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 #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 }