ESyS-Particle  4.0.1
BWallInteractionGroup.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 //       CBWallInteractionGroup functions
00015 //----------------------------------------------
00016 
00017 #include "Foundation/console.h"
00018 
00019 template<class T>
00020 CBWallInteractionGroup<T>::CBWallInteractionGroup(TML_Comm* comm):AWallInteractionGroup<T>(comm)
00021 {}
00022 
00030 template<class T>
00031 CBWallInteractionGroup<T>::CBWallInteractionGroup(TML_Comm* comm,CWall* wallp,const CBWallIGP* I)
00032   :AWallInteractionGroup<T>(comm)
00033 {
00034   console.XDebug() << "making CBWallInteractionGroup \n";
00035 
00036   m_k=I->getSpringConst();
00037   this->m_wall=wallp;
00038   m_tag=I->getTag();
00039   m_mask=I->getMask();
00040   this->m_inner_count=0;
00041 }
00042 
00043 template<class T>
00044 void CBWallInteractionGroup<T>::calcForces()
00045 {
00046 
00047   console.XDebug() << "calculating " << m_bonded_interactions.size() << " bonded wall forces\n" ;
00048   console.XDebug() << "calculating " << m_elastic_interactions.size() << " elastic wall forces\n" ;
00049 
00050   for (
00051       typename vector<CBondedWallInteraction<T> >::iterator it=m_bonded_interactions.begin();
00052       it!=m_bonded_interactions.end();
00053       it++){
00054     it->calcForces();
00055   }
00056   for(typename vector<CElasticWallInteraction<T> >::iterator it=m_elastic_interactions.begin();
00057       it!=m_elastic_interactions.end();
00058       it++){
00059     it->calcForces();
00060   }
00061 }
00062 
00063 
00070 template<class T>
00071 void CBWallInteractionGroup<T>::applyForce(const Vec3& F)
00072 {
00073   console.XDebug() << "CBWallInteractionGroup<T>::applyForce: F = " << F << "\n";
00074   // calculate local K
00075   
00076   double K=0.0;
00077   for (
00078       typename vector<CBondedWallInteraction<T> >::iterator it=m_bonded_interactions.begin();
00079       it!=m_bonded_interactions.end();
00080       it++){
00081     if(it->isInner()){
00082       K+=it->getStiffness();
00083     }
00084   }
00085   for(typename vector<CElasticWallInteraction<T> >::iterator it=m_elastic_interactions.begin();
00086       it!=m_elastic_interactions.end();
00087       it++){
00088     if(it->isInner()){
00089       K+=it->getStiffness();
00090     }
00091   }
00092   // get global K
00093   double K_global=this->m_comm->sum_all(K);
00094 
00095   int it=0;
00096   double d;
00097   Vec3 O_f=F.unit(); // direction of the applied force
00098   console.XDebug() << "CBWallInteractionGroup<T>::applyForce: unitF = " << O_f << "\n";
00099   do{
00100     // calculate local F
00101     Vec3 F_local=Vec3(0.0,0.0,0.0);
00102     // bonded itneractions
00103     for(
00104       typename vector<CBondedWallInteraction<T> >::iterator iter=m_bonded_interactions.begin();
00105       iter!=m_bonded_interactions.end();
00106       iter++
00107     ){
00108       if(iter->isInner()){
00109         Vec3 f_i=iter->getForce();
00110         F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
00111       }
00112     }
00113     // elastic interactions
00114     for(
00115       typename vector<CElasticWallInteraction<T> >::iterator iter=m_elastic_interactions.begin();
00116       iter!=m_elastic_interactions.end();
00117       iter++
00118     ){
00119       if(iter->isInner()){
00120         Vec3 f_i=iter->getForce();
00121         F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
00122       }
00123     }
00124     // get global F
00125     // by component (hack - fix later,i.e. sum_all for Vec3)
00126     double fgx=this->m_comm->sum_all(F_local.X());
00127     double fgy=this->m_comm->sum_all(F_local.Y());
00128     double fgz=this->m_comm->sum_all(F_local.Z());
00129     Vec3 F_global=Vec3(fgx,fgy,fgz);
00130 
00131     // calc necessary wall movement
00132     d=((F+F_global)*O_f)/K_global;
00133     console.XDebug()
00134       << "CBWallInteractionGroup<T>::applyForce: iteration " << it << ", d = " << fabs(d) << "\n";
00135     
00136     // move the wall
00137     console.XDebug()
00138       << "CBWallInteractionGroup<T>::applyForce: moving wall by " << d*O_f << "\n";    
00139     this->m_wall->moveBy(d*O_f);
00140     it++;
00141   } while((it<10)&&(fabs(d)>10e-6)); // check for convergence
00142   console.XDebug()
00143     << "CBWallInteractionGroup<T>::applyForce: d = " << fabs(d)
00144     << ", num iterations = " << it << "\n";
00145 }
00146 
00152 template<class T>
00153 void CBWallInteractionGroup<T>::Update(ParallelParticleArray<T>* PPA)
00154 {
00155 
00156   console.XDebug() << "CBWallInteractionGroup::Update()\n" ;
00157 
00158   // -- bonded interactions --
00159   // empty particle list first
00160   m_bonded_interactions.erase(m_bonded_interactions.begin(),m_bonded_interactions.end());
00161   this->m_inner_count=0;
00162   // build new particle list
00163   typename ParallelParticleArray<T>::ParticleListHandle plh=
00164     PPA->getParticlesAtPlane(this->m_wall->getOrigin(),this->m_wall->getNormal());
00165   for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
00166       iter!=plh->end();
00167       iter++){
00168     if(((*iter)->getTag() & m_mask )== (m_tag & m_mask)){ // if tagged->bonded ,add to bonded interaction
00169       bool iflag=PPA->isInInner((*iter)->getPos());
00170       m_bonded_interactions.push_back(CBondedWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
00171       this->m_inner_count+=(iflag ? 1 : 0);
00172     }
00173   }
00174   // -- elastic interactions --
00175   // empty particle list first
00176   m_elastic_interactions.erase(m_elastic_interactions.begin(),m_elastic_interactions.end());
00177   // build new particle list
00178   for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
00179       iter!=plh->end();
00180       iter++){
00181     if(((*iter)->getTag() & m_mask )!=(m_tag & m_mask )){ // only add to elastic if not bonded/tagged
00182       bool iflag=PPA->isInInner((*iter)->getPos());
00183       m_elastic_interactions.push_back(CElasticWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
00184       this->m_inner_count+=(iflag ? 1 : 0);
00185     }
00186   }
00187   console.XDebug() << "end CBWallInteractionGroup::Update()\n";
00188 }
00189 
00190 template<class T>
00191 ostream& operator<<(ostream& ost,const CBWallInteractionGroup<T>& IG)
00192 {
00193   ost << "CBWallInteractionGroup" << endl << flush;
00194   ost << *(IG.m_wall) << endl << flush;
00195  
00196   return ost;
00197 }