ESyS-Particle  2.3.2
BWallInteractionGroup.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 // CBWallInteractionGroup functions
15 //----------------------------------------------
16 
17 #include "Foundation/console.h"
18 
19 template<class T>
21 {}
22 
30 template<class T>
32  :AWallInteractionGroup<T>(comm)
33 {
34  console.XDebug() << "making CBWallInteractionGroup \n";
35 
36  m_k=I->getSpringConst();
37  this->m_wall=wallp;
38  m_tag=I->getTag();
39  m_mask=I->getMask();
40  this->m_inner_count=0;
41 }
42 
43 template<class T>
45 {
46 
47  console.XDebug() << "calculating " << m_bonded_interactions.size() << " bonded wall forces\n" ;
48  console.XDebug() << "calculating " << m_elastic_interactions.size() << " elastic wall forces\n" ;
49 
50  for (
51  typename vector<CBondedWallInteraction<T> >::iterator it=m_bonded_interactions.begin();
52  it!=m_bonded_interactions.end();
53  it++){
54  it->calcForces();
55  }
56  for(typename vector<CElasticWallInteraction<T> >::iterator it=m_elastic_interactions.begin();
57  it!=m_elastic_interactions.end();
58  it++){
59  it->calcForces();
60  }
61 }
62 
63 
70 template<class T>
72 {
73  console.XDebug() << "CBWallInteractionGroup<T>::applyForce: F = " << F << "\n";
74  // calculate local K
75 
76  double K=0.0;
77  for (
78  typename vector<CBondedWallInteraction<T> >::iterator it=m_bonded_interactions.begin();
79  it!=m_bonded_interactions.end();
80  it++){
81  if(it->isInner()){
82  K+=it->getStiffness();
83  }
84  }
85  for(typename vector<CElasticWallInteraction<T> >::iterator it=m_elastic_interactions.begin();
86  it!=m_elastic_interactions.end();
87  it++){
88  if(it->isInner()){
89  K+=it->getStiffness();
90  }
91  }
92  // get global K
93  double K_global=this->m_comm->sum_all(K);
94 
95  int it=0;
96  double d;
97  Vec3 O_f=F.unit(); // direction of the applied force
98  console.XDebug() << "CBWallInteractionGroup<T>::applyForce: unitF = " << O_f << "\n";
99  do{
100  // calculate local F
101  Vec3 F_local=Vec3(0.0,0.0,0.0);
102  // bonded itneractions
103  for(
104  typename vector<CBondedWallInteraction<T> >::iterator iter=m_bonded_interactions.begin();
105  iter!=m_bonded_interactions.end();
106  iter++
107  ){
108  if(iter->isInner()){
109  Vec3 f_i=iter->getForce();
110  F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
111  }
112  }
113  // elastic interactions
114  for(
115  typename vector<CElasticWallInteraction<T> >::iterator iter=m_elastic_interactions.begin();
116  iter!=m_elastic_interactions.end();
117  iter++
118  ){
119  if(iter->isInner()){
120  Vec3 f_i=iter->getForce();
121  F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
122  }
123  }
124  // get global F
125  // by component (hack - fix later,i.e. sum_all for Vec3)
126  double fgx=this->m_comm->sum_all(F_local.X());
127  double fgy=this->m_comm->sum_all(F_local.Y());
128  double fgz=this->m_comm->sum_all(F_local.Z());
129  Vec3 F_global=Vec3(fgx,fgy,fgz);
130 
131  // calc necessary wall movement
132  d=((F+F_global)*O_f)/K_global;
133  console.XDebug()
134  << "CBWallInteractionGroup<T>::applyForce: iteration " << it << ", d = " << fabs(d) << "\n";
135 
136  // move the wall
137  console.XDebug()
138  << "CBWallInteractionGroup<T>::applyForce: moving wall by " << d*O_f << "\n";
139  this->m_wall->moveBy(d*O_f);
140  it++;
141  } while((it<10)&&(fabs(d)>10e-6)); // check for convergence
142  console.XDebug()
143  << "CBWallInteractionGroup<T>::applyForce: d = " << fabs(d)
144  << ", num iterations = " << it << "\n";
145 }
146 
152 template<class T>
154 {
155 
156  console.XDebug() << "CBWallInteractionGroup::Update()\n" ;
157 
158  // -- bonded interactions --
159  // empty particle list first
160  m_bonded_interactions.erase(m_bonded_interactions.begin(),m_bonded_interactions.end());
161  this->m_inner_count=0;
162  // build new particle list
164  PPA->getParticlesAtPlane(this->m_wall->getOrigin(),this->m_wall->getNormal());
165  for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
166  iter!=plh->end();
167  iter++){
168  if(((*iter)->getTag() & m_mask )== (m_tag & m_mask)){ // if tagged->bonded ,add to bonded interaction
169  bool iflag=PPA->isInInner((*iter)->getPos());
170  m_bonded_interactions.push_back(CBondedWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
171  this->m_inner_count+=(iflag ? 1 : 0);
172  }
173  }
174  // -- elastic interactions --
175  // empty particle list first
176  m_elastic_interactions.erase(m_elastic_interactions.begin(),m_elastic_interactions.end());
177  // build new particle list
178  for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
179  iter!=plh->end();
180  iter++){
181  if(((*iter)->getTag() & m_mask )!=(m_tag & m_mask )){ // only add to elastic if not bonded/tagged
182  bool iflag=PPA->isInInner((*iter)->getPos());
183  m_elastic_interactions.push_back(CElasticWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
184  this->m_inner_count+=(iflag ? 1 : 0);
185  }
186  }
187  console.XDebug() << "end CBWallInteractionGroup::Update()\n";
188 }
189 
190 template<class T>
191 ostream& operator<<(ostream& ost,const CBWallInteractionGroup<T>& IG)
192 {
193  ost << "CBWallInteractionGroup" << endl << flush;
194  ost << *(IG.m_wall) << endl << flush;
195 
196  return ost;
197 }