ESyS-Particle  4.0.1
mesh2d_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 template<class ParticleType,class IType>
15 template<class ParticleType,class IType>
17 
25 template<class ParticleType,class IType>
27  :Mesh2D_PIS<ParticleType>(mesh_p,ppa_p),m_comm(ppa_p->getComm())
28 {
29  console.XDebug() << "Mesh2D_PIS_EB constructor\n";
30  m_param=param;
31  this->m_update_timestamp=0;
32 }
33 
43 template<class ParticleType,class IType>
44 bool Mesh2D_PIS_EB<ParticleType,IType>::isIn(const std::vector<int>& v)
45 {
46  bool res=false;
47 
48  if(v.size()<3){
49  res=false;
50  } else {
51  switch (v[2]){
52  case 1: res=m_edge_int_set.find(make_pair(v[0],v[1]))!=m_edge_int_set.end(); break;
53  case 2: res=m_edge_int_set.find(make_pair(v[0],v[1]))!=m_corner_int_set.end(); break;
54  default: console.Error() << "wrong value in argument of Mesh2D_PIS::isIn !!\n"; break;
55  }
56  }
57 
58  return res;
59 }
60 
64 template<class ParticleType,class IType>
66 {
67  console.XDebug() << "Mesh2D_PIS_EB calculating " << m_edge_interactions.size() << " edge forces and"
68  << m_corner_interactions.size() << " corner forces\n";
69 
70  // calculate forces for edge interactions
71  for(typename list<typename IType::TriIntType>::iterator ed_iter=m_edge_interactions.begin();
72  ed_iter!=m_edge_interactions.end();
73  ed_iter++){
74  ed_iter->calcForces();
75  }
76  // calculate forces for corner interactions
77  for(typename list<typename IType::CornerIntType>::iterator c_iter=m_corner_interactions.begin();
78  c_iter!=m_corner_interactions.end();
79  c_iter++){
80  c_iter->calcForces();
81  }
82 }
83 
86 template<class ParticleType,class IType>
88 {
89  console.XDebug() << "Mesh2D_PIS_EB::update on node " << m_comm.rank() << "\n";
90  bool res=false;
91 
92  // edge interactions
93  typename list<typename IType::TriIntType>::iterator iter=m_edge_interactions.begin();
94  while(iter!=m_edge_interactions.end()){
95  if(iter->broken()){
96  res=true;
97  typename list<typename IType::TriIntType>::iterator er_iter=iter;
98  iter++;
99  // remove ids from map
100  m_edge_int_set.erase(make_pair(er_iter->getTid(),er_iter->getPid()));
101  // remove interaction
102  m_edge_interactions.erase(er_iter);
103  } else {
104  iter++;
105  }
106  }
107  // corner interactions
108  typename list<typename IType::CornerIntType>::iterator c_iter=m_corner_interactions.begin();
109  while(c_iter!=m_corner_interactions.end()){
110  if(c_iter->broken()){
111  res=true;
112  typename list<typename IType::CornerIntType>::iterator cr_iter=c_iter;
113  c_iter++;
114  // remove ids from map
115  m_corner_int_set.erase(make_pair(cr_iter->getCid(),cr_iter->getPid()));
116  // remove interaction
117  m_corner_interactions.erase(cr_iter);
118  } else {
119  c_iter++;
120  }
121  }
122  console.XDebug() << "end Mesh2D_PIS_EB::update on node " << m_comm.rank() << "\n";
123  return res;
124 }
125 
126 
129 template<class ParticleType,class IType>
131 {
132  console.XDebug() << "TriMesh_PIS_EB::exchange\n";
133  for(int i=0;i<3;i++){
134  if(m_comm.get_dim(i)>1){
135  // -- up --
136  exchange_boundary(i,1);
137  // -- down --
138  exchange_boundary(i,-1);
139  }
140  }
141  console.XDebug() << "end TriMesh_PIS_EB::exchange\n";
142 }
143 
150 template<class ParticleType,class IType>
152 {
153  console.XDebug() << "Mesh2D_PIS_EB::exchange_boundary(" << dim << "," << dir << ") at node " << m_comm.rank() << "\n";
154 
155  std::set<int> bdry_ids;
156  std::vector<typename IType::TriIntType> recv_tri_buffer;
157  std::vector<typename IType::TriIntType> send_tri_buffer;
158  std::vector<typename IType::CornerIntType> recv_corner_buffer;
159  std::vector<typename IType::CornerIntType> send_corner_buffer;
160 
161  // --- setup data to send ---
162  // get boundary
163  bdry_ids = this->m_ppa->getBoundarySlabIds(dim,dir);
164  // --- edges ---
165  // for all interactions
166  for(typename std::list<typename IType::TriIntType>::iterator iter=m_edge_interactions.begin();
167  iter!=m_edge_interactions.end();
168  iter++){
169  int pid=iter->getPid();
170  if(bdry_ids.find(pid)!=bdry_ids.end()) { // if particle in interaction is in bdry -> put in buffer
171  send_tri_buffer.push_back(*iter);
172  }
173  }
174  // shift
175  m_comm.shift_cont_packed(send_tri_buffer,recv_tri_buffer,dim,dir,m_edge_exchg_tag);
176  // insert received data
177  for(typename std::vector<typename IType::TriIntType>::iterator iter=recv_tri_buffer.begin();
178  iter!=recv_tri_buffer.end();
179  iter++){
180  tryInsert(*iter);
181  }
182  // --- corners ---
183  // for all interactions
184  for(typename std::list<typename IType::CornerIntType>::iterator iter=m_corner_interactions.begin();
185  iter!=m_corner_interactions.end();
186  iter++){
187  int pid=iter->getPid();
188  if(bdry_ids.find(pid)!=bdry_ids.end()) { // if particle in interaction is in bdry -> put in buffer
189  send_corner_buffer.push_back(*iter);
190  }
191  }
192  // shift
193  m_comm.shift_cont_packed(send_corner_buffer,recv_corner_buffer,dim,dir,m_corner_exchg_tag);
194  // insert received data
195  for(typename std::vector<typename IType::CornerIntType>::iterator iter=recv_corner_buffer.begin();
196  iter!=recv_corner_buffer.end();
197  iter++){
198  tryInsert(*iter);
199  }
200 
201  console.XDebug() << "end Mesh2D_PIS_EB::exchange_boundary\n";
202 }
203 
204 template<class ParticleType,class IType>
206 {
207 }
208 
214 template<class ParticleType,class IType>
216 {
217  console.XDebug() << "Mesh2D_PIS_EB::rebuild on node " << m_comm.rank() << "\n";
219  (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be a dynamic_cast
220  // for all edge interactions
221  typename std::list<typename IType::TriIntType>::iterator ti_iter=m_edge_interactions.begin();
222  while(ti_iter!=m_edge_interactions.end()){
223  int pid=ti_iter->getPid();
224  ParticleType *part_p=t_ppa->getParticlePtrByIndex(pid);
225  if(part_p!=NULL) { // particle is available in m_ppa -> setup pointer in interaction
226  ti_iter->setPP(part_p);
227  Edge2D *ed_p = this->m_mesh->getEdgeById(ti_iter->getTid());
228  ti_iter->setTP(ed_p);
229  ti_iter++;
230  } else { // particle is not available in m_ppa -> remove interaction
231  const typename list<typename IType::TriIntType>::iterator er_iter=ti_iter;
232  ti_iter++;
233  m_edge_int_set.erase(make_pair(er_iter->getTid(),pid));
234  m_edge_interactions.erase(er_iter);
235  }
236  }
237  // and now for the corners
238  typename list<typename IType::CornerIntType>::iterator ci_iter=m_corner_interactions.begin();
239  while(ci_iter!=m_corner_interactions.end()){
240  int pid=ci_iter->getPid();
241  ParticleType *part_p=t_ppa->getParticlePtrByIndex(pid);
242  if(part_p!=NULL) { // particle is available in m_ppa -> setup pointer in interaction
243  ci_iter->setPP(part_p);
244  Corner2D *co_p = this->m_mesh->getCornerById(ci_iter->getCid());
245  ci_iter->setCP(co_p);
246  ci_iter++;
247  } else { // particle is not available in m_ppa -> remove interaction
248  const typename list<typename IType::CornerIntType>::iterator cr_iter=ci_iter;
249  ci_iter++;
250  m_corner_int_set.erase(make_pair(cr_iter->getCid(),pid));
251  m_corner_interactions.erase(cr_iter);
252  }
253  }
254 
255  console.XDebug() << "end Mesh2D_PIS_EB::rebuild on node " << m_comm.rank() << "\n";
256 }
257 
260 template<class ParticleType,class IType>
261 void Mesh2D_PIS_EB<ParticleType,IType>::tryInsert(const typename IType::TriIntType& In)
262 {
263  console.XDebug() << "Mesh2D_PIS_EB::tryInsert(const typename IType::TriIntType& In)\n";
265  (ParallelParticleArray<ParticleType>*)(this->m_ppa);
266  // check if interaction is already in
267  bool is_in=(m_edge_int_set.find(make_pair(In.getTid(),In.getPid()))!=m_edge_int_set.end());
268  ParticleType *part_p=t_ppa->getParticlePtrByIndex(In.getPid());
269  if((!is_in) && (part_p!=NULL)){
270  m_edge_interactions.push_back(In);
271  m_edge_int_set.insert(make_pair(In.getTid(),In.getPid()));
272  }
273 }
274 
277 template<class ParticleType,class IType>
278 void Mesh2D_PIS_EB<ParticleType,IType>::tryInsert(const typename IType::CornerIntType& In)
279 {
280  console.XDebug() << "Mesh2D_PIS_EB::tryInsert(const typename IType::CornerIntType& In)\n";
282  (ParallelParticleArray<ParticleType>*)(this->m_ppa);
283  // check if interaction is already in
284  bool is_in=(m_corner_int_set.find(make_pair(In.getCid(),In.getPid()))!=m_corner_int_set.end());
285  ParticleType *part_p=t_ppa->getParticlePtrByIndex(In.getPid());
286  if((!is_in) && (part_p!=NULL)){
287  m_corner_interactions.push_back(In);
288  m_corner_int_set.insert(make_pair(In.getCid(),In.getPid()));
289  }
290 }
291 
303 template<class ParticleType,class IType>
304 void Mesh2D_PIS_EB<ParticleType,IType>::tryInsert(const vector<int>& pids)
305 {
306  console.XDebug() << "Mesh2D_PIS_EB::(const vector<int>& pids)\n";
308  (ParallelParticleArray<ParticleType>*)(this->m_ppa); // should be dynamic_cast
309 
310  if(pids.size()<3){
311  bool is_in=isIn(pids); // interaction already in
312  ParticleType *part_p=t_ppa->getParticlePtrByIndex(pids[1]);
313  if((!is_in) && (part_p!=NULL)){
314  switch (pids[2]){
315  case 1: { // edge
316  Edge2D *edge_p = this->m_mesh->getEdgeById(pids[0]);
317  if(edge_p!=NULL){
318  m_edge_interactions.push_back(
319  typename IType::TriIntType(
320  part_p,
321  edge_p,
322  m_param,
323  this->m_ppa->isInInner(part_p->getPos())
324  )
325  );
326  m_edge_int_set.insert(make_pair(pids[0],pids[1]));
327  } else {
328  console.Error() << "ERROR: Wrong edge id " << pids[0] << " in Mesh2D_PIS_EB::tryInsert\n";
329  }
330  } break;
331  case 2: {
332  Corner2D *corner_p = this->m_mesh->getCornerById(pids[0]);
333  if(corner_p!=NULL){
334  m_corner_interactions.push_back(
335  typename IType::CornerIntType(
336  part_p,
337  corner_p,
338  m_param,
339  this->m_ppa->isInInner(part_p->getPos())
340  )
341  );
342  m_corner_int_set.insert(make_pair(pids[0],pids[1]));
343  } else {
344  console.Error() << "ERROR: Wrong corner id " << pids[0] << " in Mesh2D_PIS_EB::tryInsert\n";
345  }
346  } break;
347  default : console.Error() << "Error: pids[2]= " << pids[2] << "\n"; // Can't happen !!
348  }
349  }
350  }
351 }
352 
353 
356 template<class ParticleType,class IType>
358 {
359  console.XDebug() << "Mesh2D_PIS_EB::buildFromPPATagged(" << tag << "," << mask << ") not implemented \n";
360 
361  console.XDebug() << "end Mesh2D_PIS_EB::buildFromPPATagged()";
362 }
363 
369 template<class ParticleType,class IType>
371 {
372  console.XDebug() << "Mesh2D_PIS_EB::buildFromPPAByGap(" << gmax << ")\n";
373  set<int> id_set;
374 
375  // for all edges
376  for(
377  Mesh2D::edge_iterator ed_iter = this->m_mesh->edges_begin();
378  ed_iter != this->m_mesh->edges_end();
379  ed_iter++
380  ){
381  // get particles near edge
382  typename ParallelParticleArray<ParticleType>::ParticleListHandle plh=
383  ((ParallelParticleArray<ParticleType>*)this->m_ppa)->getParticlesNearEdge(&(*ed_iter));
384  console.Debug() << "edge " << ed_iter->getID() << " nr. of particles : " << plh->size() << "\n";
385  // for all particles found
386  for(
387  typename ParallelParticleArray<ParticleType>::ParticleListIterator p_iter=plh->begin();
388  p_iter!=plh->end();
389  p_iter++){
390  // if valid interaction
391  console.Debug() << "interaction : " << ed_iter->getID() << " " << (*p_iter)->getID() << "\n";
392  if(id_set.find((*p_iter)->getID())==id_set.end()){
393  pair<bool,double> dist=ed_iter->dist((*p_iter)->getPos());
394  console.Debug() << "is valid: " << dist.first << " dist : " << dist.second << "\n";
395  if(dist.first){
396  // check ga
397  double gap=fabs(dist.second-(*p_iter)->getRad());
398  console.Debug() << "radius: " << (*p_iter)->getRad() << " gap : " << gap << "\n";
399  // if small enough, add
400  if(gap<gmax){
401  console.Debug() << "Inserting " << (*p_iter)->getID() << " !!\n";
402  bool in_flag = this->m_ppa->isInInner((*p_iter)->getPos());
403  m_edge_interactions.push_back(typename IType::TriIntType((*p_iter),&(*ed_iter),m_param,in_flag));
404  m_edge_int_set.insert(make_pair(ed_iter->getID(),(*p_iter)->getID()));
405  id_set.insert((*p_iter)->getID());
406  }
407  }
408  }
409  }
410  }
411 }
412 
413 template <class ParticleType,class IType>
415 {
416  return
417  InteractionIterator(
418  m_edge_interactions.begin(),
419  m_edge_interactions.end(),
420  this->m_ppa
421  );
422 }
423 
424 template <class ParticleType,class IType>
426 {
427  const std::string delim = "\n";
428  typedef typename IType::TriIntType::CheckPointable CheckPointable;
429 
430  InteractionIterator it = getInnerInteractionIterator();
431  ost << IType::getType() << delim;
432  ost << it.getNumRemaining();
433  while (it.hasNext()){
434  ost << delim;
435  CheckPointable(it.next()).saveCheckPointData(ost);
436  }
437 }
438 
439 template <class ParticleType,class IType>
441 {
442  console.Critical() << "Mesh2D_PIS_EB<ParticleType,IType>::loadCheckPointData NOT IMPLEMENTED\n";
443 }
444 
445 // --- Iterator members ---
446 
447 template <class ParticleType,class IType>
449  : m_numRemaining(0),
450  m_it(end),
451  m_end(end),
452  m_ppa(ppa)
453 {
454  m_numRemaining = 0;
455  for (Iterator it = begin; it != end; it++) {
456  if (isInner(it)) {
457  m_numRemaining++;
458  }
459  }
460  m_it = begin;
461  m_end = end;
462 }
463 
464 
465 template <class ParticleType,class IType>
467 {
468  return (m_numRemaining > 0);
469 }
470 
471 template <class ParticleType,class IType>
473 {
474  while (!isInner(m_it)) {
475  m_it++;
476  }
477  Interaction &i = *m_it;
478  m_it++;
479  m_numRemaining--;
480  return i;
481 }
482 
483 template <class ParticleType,class IType>
485 {
486  return m_numRemaining;
487 }
488 
489 template <class ParticleType,class IType>
491 {
492  return m_ppa->isInInner(it->getPos());
493 }