ESyS-Particle  4.0.1
ViscWallIG.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 //----------------------------------------------
14 // CViscWallIG functions
15 //----------------------------------------------
16 
17 #include "Foundation/console.h"
18 
19 template<class T>
21 {
22 }
23 
31 template<class T>
33  :AWallInteractionGroup<T>(comm)
34 {
35  console.XDebug() << "making CViscWallIG pos \n";
36 
37  m_k=I->getSpringConst();
38  this->m_wall=wallp;
39  m_tag=I->getTag();
40  m_nu=I->getNu();
41  this->m_inner_count=0;
42 }
43 
44 template<class T>
46 {
47 
48  console.XDebug() << "calculating " << m_visc_interactions.size() << " viscous wall forces\n" ;
49  console.XDebug() << "calculating " << m_elastic_interactions.size() << " elastic wall forces\n" ;
50 
51  for(
52  typename vector<CViscWallInteraction<T> >::iterator it=m_visc_interactions.begin();
53  it!=m_visc_interactions.end();
54  it++
55  ){
56  it->calcForces();
57  }
58  for(
59  typename vector<CElasticWallInteraction<T> >::iterator it=m_elastic_interactions.begin();
60  it!=m_elastic_interactions.end();
61  it++
62  ){
63  it->calcForces();
64  }
65 }
66 
72 template<class T>
74 {
75  this->m_wall->setVel(V);
76 }
77 
84 template<class T>
86 {
87  // calculate local K
88  double K=this->m_inner_count*m_k;
89  // get global K
90  double K_global=this->m_comm->sum_all(K);
91 
92  int it=0;
93  double d;
94  Vec3 O_f=F.unit(); // direction of the applied force
95  do{
96  // calculate local F
97  Vec3 F_local=Vec3(0.0,0.0,0.0);
98  // viscous interactions
99  for (
100  typename vector<CViscWallInteraction<T> >::iterator iter=m_visc_interactions.begin();
101  iter != m_visc_interactions.end();
102  iter++
103  ){
104  if(iter->isInner()){
105  Vec3 f_i=iter->getForce();
106  F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
107  }
108  }
109  // elastic interactions
110  for (
111  typename vector<CElasticWallInteraction<T> >::iterator iter=m_elastic_interactions.begin();
112  iter != m_elastic_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)/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 
140 template<class T>
142 {
143 
144  console.XDebug() << "CViscWallIG::Update()\n" ;
145 
146  // -- bonded interactions --
147  // empty particle list first
148  m_visc_interactions.erase(m_visc_interactions.begin(),m_visc_interactions.end());
149  this->m_inner_count=0;
150  // build new particle list
151  typename ParallelParticleArray<T>::ParticleListHandle plh=
152  PPA->getParticlesAtPlane(this->m_wall->getOrigin(),this->m_wall->getNormal());
153  for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
154  iter!=plh->end();
155  iter++){
156  if((*iter)->getTag()==m_tag){// if tagged -> apply viscous drag, i.e. add to viscous interaction
157  bool iflag=PPA->isInInner((*iter)->getPos());
158  m_visc_interactions.push_back(CViscWallInteraction<T>(*iter,this->m_wall,m_nu,iflag));
159  this->m_inner_count+=(iflag ? 1 : 0);
160  }
161  }
162  // -- elastic interactions --
163  // empty particle list first
164  m_elastic_interactions.erase(m_elastic_interactions.begin(),m_elastic_interactions.end());
165  // build new particle list
166  for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
167  iter!=plh->end();
168  iter++){ // always add to elastic group, even if in viscous
169  bool iflag=PPA->isInInner((*iter)->getPos());
170  m_elastic_interactions.push_back(CElasticWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
171  this->m_inner_count+=(iflag ? 1 : 0);
172  }
173  console.XDebug() << "end CViscWallIG::Update()\n";
174 }
175 
176 template<class T>
177 ostream& operator<<(ostream& ost,const CViscWallIG<T>& IG)
178 {
179  ost << "CViscWallIG" << endl << flush;
180  ost << *(IG.m_wall) << endl << flush;
181 
182  return ost;
183 }