ESyS-Particle
4.0.1
|
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 }