ESyS-Particle  4.0.1
ViscWallIG.hpp
00001 
00002 //                                                         //
00003 // Copyright (c) 2003-2011 by The University of Queensland //
00004 // Earth Systems Science Computational Centre (ESSCC)      //
00005 // http://www.uq.edu.au/esscc                              //
00006 //                                                         //
00007 // Primary Business: Brisbane, Queensland, Australia       //
00008 // Licensed under the Open Software License version 3.0    //
00009 // http://www.opensource.org/licenses/osl-3.0.php          //
00010 //                                                         //
00012 
00013 //----------------------------------------------
00014 //       CViscWallIG functions
00015 //----------------------------------------------
00016 
00017 #include "Foundation/console.h"
00018 
00019 template<class T>
00020 CViscWallIG<T>::CViscWallIG(TML_Comm* comm):AWallInteractionGroup<T>(comm)
00021 {
00022 }
00023 
00031 template<class T>
00032 CViscWallIG<T>::CViscWallIG(TML_Comm* comm,CWall* wallp,const CVWallIGP* I)
00033   :AWallInteractionGroup<T>(comm)
00034 {
00035   console.XDebug() << "making CViscWallIG pos \n";
00036 
00037   m_k=I->getSpringConst();
00038   this->m_wall=wallp;
00039   m_tag=I->getTag();
00040   m_nu=I->getNu();
00041   this->m_inner_count=0;
00042 }
00043 
00044 template<class T>
00045 void CViscWallIG<T>::calcForces()
00046 {
00047 
00048   console.XDebug() << "calculating " << m_visc_interactions.size() << " viscous wall forces\n" ;
00049   console.XDebug() << "calculating " << m_elastic_interactions.size() << " elastic wall forces\n" ;
00050 
00051   for(
00052     typename vector<CViscWallInteraction<T> >::iterator it=m_visc_interactions.begin();
00053     it!=m_visc_interactions.end();
00054     it++
00055   ){
00056     it->calcForces();
00057   }
00058   for(
00059     typename vector<CElasticWallInteraction<T> >::iterator it=m_elastic_interactions.begin();
00060     it!=m_elastic_interactions.end();
00061     it++
00062   ){
00063     it->calcForces();
00064   }
00065 }
00066 
00072 template<class T>
00073 void CViscWallIG<T>::setVelocity(const Vec3& V)
00074 {
00075   this->m_wall->setVel(V);
00076 }
00077 
00084 template<class T>
00085 void CViscWallIG<T>::applyForce(const Vec3& F)
00086 {
00087   // calculate local K
00088   double K=this->m_inner_count*m_k;
00089   // get global K
00090   double K_global=this->m_comm->sum_all(K);
00091 
00092   int it=0;
00093   double d;
00094   Vec3 O_f=F.unit(); // direction of the applied force
00095   do{
00096     // calculate local F
00097     Vec3 F_local=Vec3(0.0,0.0,0.0);
00098     // viscous interactions
00099     for (
00100       typename vector<CViscWallInteraction<T> >::iterator iter=m_visc_interactions.begin();
00101         iter != m_visc_interactions.end();
00102         iter++
00103       ){
00104       if(iter->isInner()){
00105         Vec3 f_i=iter->getForce();
00106         F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
00107       }
00108     }
00109     // elastic interactions
00110     for (
00111       typename vector<CElasticWallInteraction<T> >::iterator iter=m_elastic_interactions.begin();
00112       iter != m_elastic_interactions.end();
00113       iter++
00114     ){
00115       if(iter->isInner()){
00116         Vec3 f_i=iter->getForce();
00117         F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
00118       }
00119     }
00120     // get global F
00121     // by component (hack - fix later,i.e. sum_all for Vec3)
00122     double fgx=this->m_comm->sum_all(F_local.X());
00123     double fgy=this->m_comm->sum_all(F_local.Y());
00124     double fgz=this->m_comm->sum_all(F_local.Z());
00125     Vec3 F_global=Vec3(fgx,fgy,fgz);
00126 
00127     // calc necessary wall movement
00128     d=((F+F_global)*O_f)/K_global;
00129     // move the wall
00130     this->m_wall->moveBy(d*O_f);
00131     it++;
00132   } while((it<10)&&(fabs(d)>10e-6)); // check for convergence
00133 }
00134 
00140 template<class T>
00141 void CViscWallIG<T>::Update(ParallelParticleArray<T>* PPA)
00142 {
00143 
00144   console.XDebug() << "CViscWallIG::Update()\n" ;
00145 
00146   // -- bonded interactions --
00147   // empty particle list first
00148   m_visc_interactions.erase(m_visc_interactions.begin(),m_visc_interactions.end());
00149   this->m_inner_count=0;
00150   // build new particle list
00151   typename ParallelParticleArray<T>::ParticleListHandle plh=
00152     PPA->getParticlesAtPlane(this->m_wall->getOrigin(),this->m_wall->getNormal());
00153   for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
00154       iter!=plh->end();
00155       iter++){
00156     if((*iter)->getTag()==m_tag){// if tagged -> apply viscous drag, i.e. add to viscous interaction
00157       bool iflag=PPA->isInInner((*iter)->getPos());
00158       m_visc_interactions.push_back(CViscWallInteraction<T>(*iter,this->m_wall,m_nu,iflag));
00159       this->m_inner_count+=(iflag ? 1 : 0);
00160     }
00161   }
00162   // -- elastic interactions --
00163   // empty particle list first
00164   m_elastic_interactions.erase(m_elastic_interactions.begin(),m_elastic_interactions.end());
00165   // build new particle list
00166   for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
00167       iter!=plh->end();
00168       iter++){ // always add to elastic group, even if in viscous
00169     bool iflag=PPA->isInInner((*iter)->getPos());
00170     m_elastic_interactions.push_back(CElasticWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
00171     this->m_inner_count+=(iflag ? 1 : 0);
00172   }
00173   console.XDebug() << "end CViscWallIG::Update()\n";
00174 }
00175 
00176 template<class T>
00177 ostream& operator<<(ostream& ost,const CViscWallIG<T>& IG)
00178 {
00179   ost << "CViscWallIG" << endl << flush;
00180   ost << *(IG.m_wall) << endl << flush;
00181  
00182   return ost;
00183 }