ESyS-Particle  2.3.2
brokenEWallInteractionGroup.hpp
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
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 //----------------------------------------------
14 // CEWallInteractionGroup functions
15 //----------------------------------------------
16 
17 #include "Foundation/console.h"
18 #include <iostream>
19 
20 template<class T>
22 {}
23 
31 template<class T>
33  :AWallInteractionGroup<T>(comm)
34 {
35  console.XDebug() << "making CEWallInteractionGroup \n";
36 
37  m_k=I->getSpringConst();
38  this->m_wall=wallp;
39 }
40 
41 template<class T>
43 {
44 
45  console.XDebug() << "calculating " << m_interactions.size() << " elastic wall forces\n" ;
46 
47  for(
48  typename vector<CElasticWallInteraction<T> >::iterator it=m_interactions.begin();
49  it != m_interactions.end();
50  it++
51  ){
52  it->calcForces();
53  }
54 }
55 
56 template<class T>
58 {
59 
60  console.XDebug() << "CEWallInteractionGroup::Update()\n" ;
61 
62  console.XDebug()
63  << "CEWallInteractionGroup::Update: wall origin = " << this->m_wall->getOrigin()
64  << ", wall normal = " << this->m_wall->getNormal() << "\n" ;
65 
66  double k_local=0.0;
67  // empty particle list first
68  m_interactions.erase(m_interactions.begin(),m_interactions.end());
69  this->m_inner_count=0;
70  // build new particle list
72  PPA->getParticlesAtPlane(this->m_wall->getOrigin(),this->m_wall->getNormal());
73  for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
74  iter!=plh->end();
75  iter++){
76  bool iflag=PPA->isInInner((*iter)->getPos());
77  m_interactions.push_back(CElasticWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
78  this->m_inner_count+=(iflag ? 1 : 0);
79  if(iflag){
81  k_local+=m_k*((*iter)->getRad()); // update local K
82  } else {
83  k_local+=m_k;
84  }
85  }
86  }
87  // get global K
88  m_k_global=this->m_comm->sum_all(k_local);
89 
90  console.XDebug() << "end CEWallInteractionGroup::Update()\n";
91 }
92 
93 
101 template<class T>
103 {
104  int it=0;
105  double d;
106  Vec3 O_f=F.unit(); // direction of the applied force
107  do{
108  // calculate local F
109  Vec3 F_local=Vec3(0.0,0.0,0.0);
110  for (
111  typename vector<CElasticWallInteraction<T> >::iterator iter=m_interactions.begin();
112  iter!=m_interactions.end();
113  iter++
114  ){
115  if(iter->isInner()){
116  Vec3 f_i=iter->getForce();
117  F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
118  }
119  }
120  // get global F
121  // by component (hack - fix later,i.e. sum_all for Vec3)
122  double fgx=this->m_comm->sum_all(F_local.X());
123  double fgy=this->m_comm->sum_all(F_local.Y());
124  double fgz=this->m_comm->sum_all(F_local.Z());
125  Vec3 F_global=Vec3(fgx,fgy,fgz);
126 
127  // calc necessary wall movement
128  d=((F+F_global)*O_f)/m_k_global;
129  // move the wall
130  this->m_wall->moveBy(d*O_f);
131  it++;
132  } while((it<10)&&(fabs(d)>10e-6)); // check for convergence
133 }
134 
135 
136 template<class T>
137 ostream& operator<<(ostream& ost,const CEWallInteractionGroup<T>& IG)
138 {
139  ost << "CEWallInteractionGroup" << endl << flush;
140  ost << *(IG.m_wall) << endl << flush;
141 
142  return ost;
143 }