ESyS-Particle  4.0.1
EWallInteractionGroup.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 // 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  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
71  typename ParallelParticleArray<T>::ParticleListHandle plh=
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  }
80 
81  console.XDebug() << "end CEWallInteractionGroup::Update()\n";
82 }
83 
84 
92 template<class T>
94 {
95 
96  int it=0;
97  double d; // distance to move
98  double df; // force difference
99  double ef; // relative force error (df/|F|)
100  Vec3 O_f=F.unit(); // direction of the applied force
101  do{
102  //std::cerr << "iteration: " << it << std::endl;
103  // calculate local stiffness
104  k_local=0.0;
105  for(typename vector<CElasticWallInteraction<T> >::iterator iter=m_interactions.begin();
106  iter!=m_interactions.end();
107  iter++){
108  k_local+=iter->getStiffness();
109  }
110  //std::cerr << "local stiffness: " << k_local << std::endl;
111  // get global K
112  m_k_global=this->m_comm->sum_all(k_local);
113  //std::cerr << "global stiffness: " << m_k_global << std::endl;
114  if(m_k_global>0){
115  // calculate local F
116  Vec3 F_local=Vec3(0.0,0.0,0.0);
117  for (
118  typename vector<CElasticWallInteraction<T> >::iterator iter=m_interactions.begin();
119  iter!=m_interactions.end();
120  iter++
121  ){
122  if(iter->isInner()){
123  Vec3 f_i=iter->getForce();
124  F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
125  }
126  }
127  //std::cerr << "local Force: " << F_local << std::endl;
128  // get global F
129  // by component (hack - fix later,i.e. sum_all for Vec3)
130  double fgx=this->m_comm->sum_all(F_local.X());
131  double fgy=this->m_comm->sum_all(F_local.Y());
132  double fgz=this->m_comm->sum_all(F_local.Z());
133  Vec3 F_global=Vec3(fgx,fgy,fgz);
134  //std::cerr << "global Force: " << F_global << std::endl;
135 
136  // calc necessary wall movement
137  df=(F+F_global)*O_f;
138  d=df/m_k_global;
139  ef=df/F.norm();
140  //std::cerr << "d,df,ef: " << d << " , " << df << " , " << ef << std::endl;
141  // move the wall
142  //std::cerr << " old wall pos: " << this->m_wall->getPos() << std::endl;
143  this->m_wall->moveBy(d*O_f);
144  //std::cerr << "new wall pos: " << this->m_wall->getPos() << std::endl;
145  it++;
146  } else {
147  d=1e-5;
148  ef=1;
149  //std::cerr << "no contact" << std::endl;
150  //std::cerr << "old wall pos: " << this->m_wall->getPos() << std::endl;
151  this->m_wall->moveBy(d*O_f);
152  //std::cerr << "new wall pos: " << this->m_wall->getPos() << std::endl;
153  }
154  } while((it<50)&&(ef>1e-3)); // check for convergence
155 }
156 
157 
158 template<class T>
159 ostream& operator<<(ostream& ost,const CEWallInteractionGroup<T>& IG)
160 {
161  ost << "CEWallInteractionGroup" << endl << flush;
162  ost << *(IG.m_wall) << endl << flush;
163 
164  return ost;
165 }