ESyS-Particle  4.0.1
trimesh_pis_eb.hpp
1 
2 // //
3 // Copyright (c) 2003-2011 by The University of Queensland //
4 // Earth Systems Science Computational Centre (ESSCC) //
5 // http://www.uq.edu.au/esscc //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 #include "Foundation/console.h"
14 
15 template<class ParticleType,class IType>
17 
25 template<class ParticleType,class IType>
27  : TriMesh_PIS<ParticleType>(mesh_p,ppa_p),m_comm(ppa_p->getComm())
28 {
29  console.XDebug() << "TriMesh_PIS_EB constructor\n";
30  m_param=param;
31  this->m_update_timestamp = 0;
32 }
33 
34 template <class ParticleType,class IType>
36 {
37 }
38 
48 template <class ParticleType,class IType>
49 bool TriMesh_PIS_EB<ParticleType,IType>::isIn(const std::vector<int>& v)
50 {
51  bool res=false;
52 
53  if(v.size()<3){
54  res=false;
55  } else {
56  switch (v[2]){
57  case 0: res=m_tri_int_set.find(make_pair(v[0],v[1]))!=m_tri_int_set.end(); break;
58  default: console.Error() << "wrong value in argument of TriMesh_PIS::isIn !!\n"; break;
59  }
60  }
61 
62  return res;
63 }
64 
68 template<class ParticleType,class IType>
70 {
71  console.XDebug() << "TriMesh_PIS_EB calculating " << m_triangle_interactions.size() << " triangle forces\n";
72 
73  // calculate forces for triangle interactions
74  for(typename list<typename IType::TriIntType>::iterator tri_iter=m_triangle_interactions.begin();
75  tri_iter!=m_triangle_interactions.end();
76  tri_iter++){
77  tri_iter->calcForces();
78  }
79 }
80 
83 template<class ParticleType,class IType>
85 {
86  console.XDebug() << "TriMesh_PIS_EB::update on node " << m_comm.rank() << "\n";
87  bool res=false;
88 
89  typename list<typename IType::TriIntType>::iterator iter=m_triangle_interactions.begin();
90  while(iter!=m_triangle_interactions.end()){
91  if(iter->broken()){
92  res=true;
93  typename list<typename IType::TriIntType>::iterator er_iter=iter;
94  iter++;
95  // remove ids from map
96  m_tri_int_set.erase(make_pair(er_iter->getTid(),er_iter->getPid()));
97  // remove interaction
98  m_triangle_interactions.erase(er_iter);
99  } else {
100  iter++;
101  }
102  }
103 
104  console.XDebug() << "end TriMesh_PIS_EB::update on node " << m_comm.rank() << "\n";
105  return res;
106 }
107 
114 template<class ParticleType,class IType>
116 {
117  console.XDebug() << "TriMesh_PIS_EB::exchange_boundary(" << dim << "," << dir << ") at node " << m_comm.rank() << "\n";
118 
119  std::set<int> bdry_ids;
120  std::vector<typename IType::TriIntType> recv_tri_buffer;
121  std::vector<typename IType::TriIntType> send_tri_buffer;
122 
123  // --- setup data to send ---
124  // get boundary
125  bdry_ids = this->m_ppa->getBoundarySlabIds(dim,dir);
126  // for all interactions
127  for(typename list<typename IType::TriIntType>::iterator iter=m_triangle_interactions.begin();
128  iter!=m_triangle_interactions.end();
129  iter++){
130  int pid=iter->getPid();
131  if(bdry_ids.find(pid)!=bdry_ids.end()) { // if particle in interaction is in bdry -> put in buffer
132  send_tri_buffer.push_back(*iter);
133  }
134  }
135  // ---- shift ----
136  m_comm.shift_cont_packed(send_tri_buffer,recv_tri_buffer,dim,dir,m_exchg_tag);
137  // insert received data
138  for(typename std::vector<typename IType::TriIntType>::iterator iter=recv_tri_buffer.begin();
139  iter!=recv_tri_buffer.end();
140  iter++){
141  tryInsert(*iter);
142  }
143 
144  console.XDebug() << "end TriMesh_PIS_EB::exchange_boundary\n";
145 }
146 
149 template<class ParticleType,class IType>
151 {
152  console.XDebug() << "TriMesh_PIS_EB::exchange\n";
153  for(int i=0;i<3;i++){
154  if(m_comm.get_dim(i)>1){
155  // -- up --
156  exchange_boundary(i,1);
157  // -- down --
158  exchange_boundary(i,-1);
159  }
160  }
161  console.XDebug() << "end TriMesh_PIS_EB::exchange\n";
162 }
163 
169 template<class ParticleType,class IType>
171 {
172  console.XDebug() << "TriMesh_PIS_EB::rebuild on node " << m_comm.rank() << "\n";
174  (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be a dynamic_cast
175  // for all triangle interactions
176  typename list<typename IType::TriIntType>::iterator ti_iter=m_triangle_interactions.begin();
177  while(ti_iter!=m_triangle_interactions.end()){
178  int pid=ti_iter->getPid();
179  ParticleType *part_p=t_ppa->getParticlePtrByIndex(pid);
180  if(part_p!=NULL) { // particle is available in m_ppa -> setup pointer in interaction
181  ti_iter->setPP(part_p);
182  Triangle *tri_p = this->m_mesh->getTriangleById(ti_iter->getTid());
183  ti_iter->setTP(tri_p);
184  ti_iter++;
185  } else { // particle is not available in m_ppa -> remove interaction
186  const typename list<typename IType::TriIntType>::iterator er_iter=ti_iter;
187  ti_iter++;
188  m_tri_int_set.erase(make_pair(er_iter->getTid(),pid));
189  m_triangle_interactions.erase(er_iter);
190  }
191  }
192 
193  console.XDebug() << "end TriMesh_PIS_EB::rebuild on node " << m_comm.rank() << "\n";
194 }
195 
207 template<class ParticleType,class IType>
208 void TriMesh_PIS_EB<ParticleType,IType>::tryInsert(const std::vector<int>& pids)
209 {
211  (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be dynamic_cast
212 
213  if(pids.size()<3){
214  bool is_in=isIn(pids); // interaction already in
215  ParticleType *part_p=t_ppa->getParticlePtrByIndex(pids[1]);
216  if((!is_in) && (part_p!=NULL)){
217  switch (pids[2]){
218  case 0: {
219  Triangle *tri_p = this->m_mesh->getTriangleById(pids[0]);
220  if (tri_p!=NULL){
221  m_triangle_interactions.push_back(
222  typename IType::TriIntType(
223  part_p,
224  tri_p,
225  m_param,
226  this->m_ppa->isInInner(part_p->getPos())
227  )
228  );
229  m_tri_int_set.insert(make_pair(pids[0],pids[1]));
230  }
231  } break;
232  default : {
233  console.Error()
234  << "Error: pids[2]= " << pids[2]
235  << "\n"; // Can't happen !!
236  }
237  }
238  }
239  }
240 }
241 
244 template<class ParticleType,class IType>
245 void TriMesh_PIS_EB<ParticleType,IType>::tryInsert(const typename IType::TriIntType& In)
246 {
248  (ParallelParticleArray<ParticleType>*)(this->m_ppa);
249  // check if interaction is already in
250  bool is_in=(m_tri_int_set.find(make_pair(In.getTid(),In.getPid()))!=m_tri_int_set.end());
251  ParticleType *part_p=t_ppa->getParticlePtrByIndex(In.getPid());
252  if((!is_in) && (part_p!=NULL)){
253  m_triangle_interactions.push_back(In);
254  m_tri_int_set.insert(make_pair(In.getTid(),In.getPid()));
255  }
256 }
257 
260 template<class ParticleType,class IType>
262 {
263  console.XDebug() << "TriMesh_PIS_EB::buildFromPPATagged(" << tag << "," << mask << ")\n";
264 
265  console.XDebug() << "end TriMesh_PIS_EB::buildFromPPATagged()";
266 }
267 
270 template<class ParticleType,class IType>
272 {
273  console.XDebug() << "TriMesh_PIS_EB::buildFromPPAByGap(" << gmax << ")\n";
274  set<int> id_set;
275 
276  // for all triangles
277  for (
278  TriMesh::triangle_iterator tri_iter = this->m_mesh->triangles_begin();
279  tri_iter != this->m_mesh->triangles_end();
280  tri_iter++
281  ){
282  // get particles near triangle
283  typename ParallelParticleArray<ParticleType>::ParticleListHandle plh=
284  ((ParallelParticleArray<ParticleType>*)this->m_ppa)->getParticlesNearTriangle(*tri_iter);
285  console.XDebug() << "triangle " << tri_iter->getID() << " nr. of particles : " << plh->size() << "\n";
286  // for all particles found
287  for(typename ParallelParticleArray<ParticleType>::ParticleListIterator p_iter=plh->begin();
288  p_iter!=plh->end();
289  p_iter++){
290  // if valid interaction
291  console.XDebug() << "interaction : " << tri_iter->getID() << " " << (*p_iter)->getID() << "\n";
292  if(id_set.find((*p_iter)->getID())==id_set.end()){
293  pair<bool,double> dist=tri_iter->dist((*p_iter)->getPos());
294  console.XDebug() << "is valid: " << dist.first << " dist : " << dist.second << "\n";
295  if(dist.first){
296  // check gap
297  double gap=fabs(dist.second-(*p_iter)->getRad());
298  console.XDebug() << "radius: " << (*p_iter)->getRad() << " gap : " << gap << "\n";
299  // if small enough, add
300  if(gap<gmax){
301  console.XDebug() << "Insert !!\n";
302  bool in_flag = this->m_ppa->isInInner((*p_iter)->getPos());
303  m_triangle_interactions.push_back(typename IType::TriIntType((*p_iter),&(*tri_iter),m_param,in_flag));
304  m_tri_int_set.insert(make_pair(tri_iter->getID(),(*p_iter)->getID()));
305  id_set.insert((*p_iter)->getID());
306  }
307  }
308  }
309  }
310  }
311  console.XDebug() << "end TriMesh_PIS_EB::buildFromPPAByGap()";
312 }
313 
319 template<class ParticleType,class IType>
321 {
322  const std::string delim = "\n";
323  typedef typename IType::TriIntType::CheckPointable CheckPointable;
324 
325  // stage 1: count how many interactions with "inner" particles we have
326  int icount=0;
327  for(typename list<typename IType::TriIntType>::iterator it=m_triangle_interactions.begin();
328  it!=m_triangle_interactions.end();
329  it++){
330  if(it->isInner()) icount++;
331  }
332 
333  // stage 2: write data
334  oStream << IType::getType() << delim;
335  oStream << icount << delim;
336  for(typename list<typename IType::TriIntType>::iterator it=m_triangle_interactions.begin();
337  it!=m_triangle_interactions.end();
338  it++){
339  if(it->isInner()) CheckPointable(*it).saveCheckPointData(oStream);
340  oStream << delim;
341  }
342 }