ESyS-Particle  4.0.1
EWallInteractionGroup.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 //       CEWallInteractionGroup functions
00015 //----------------------------------------------
00016 
00017 #include "Foundation/console.h"
00018 #include <iostream>
00019 
00020 template<class T>
00021 CEWallInteractionGroup<T>::CEWallInteractionGroup(TML_Comm* comm):AWallInteractionGroup<T>(comm)
00022 {}
00023 
00031 template<class T>
00032 CEWallInteractionGroup<T>::CEWallInteractionGroup(TML_Comm* comm,CWall* wallp,const CEWallIGP* I)
00033   :AWallInteractionGroup<T>(comm)
00034 {
00035   console.XDebug() << "making CEWallInteractionGroup \n";
00036 
00037   m_k=I->getSpringConst();
00038   this->m_wall=wallp;
00039 }
00040 
00041 template<class T>
00042 void CEWallInteractionGroup<T>::calcForces()
00043 {
00044 
00045   console.XDebug() << "calculating " << m_interactions.size() << " elastic wall forces\n" ;
00046 
00047   for(
00048     typename vector<CElasticWallInteraction<T> >::iterator it=m_interactions.begin();
00049     it != m_interactions.end();
00050     it++
00051   ){
00052     it->calcForces();
00053   }
00054 }
00055 
00056 template<class T>
00057 void CEWallInteractionGroup<T>::Update(ParallelParticleArray<T>* PPA)
00058 {
00059 
00060   console.XDebug() << "CEWallInteractionGroup::Update()\n" ;
00061 
00062   console.XDebug()
00063     << "CEWallInteractionGroup::Update: wall origin = " << this->m_wall->getOrigin()
00064     << ", wall normal = " << this->m_wall->getNormal() << "\n" ;
00065 
00066   k_local=0.0;
00067   // empty particle list first
00068   m_interactions.erase(m_interactions.begin(),m_interactions.end());
00069   this->m_inner_count=0;
00070   // build new particle list
00071   typename ParallelParticleArray<T>::ParticleListHandle plh=
00072     PPA->getParticlesAtPlane(this->m_wall->getOrigin(),this->m_wall->getNormal());
00073   for(typename ParallelParticleArray<T>::ParticleListIterator iter=plh->begin();
00074       iter!=plh->end();
00075       iter++){
00076     bool iflag=PPA->isInInner((*iter)->getPos());
00077     m_interactions.push_back(CElasticWallInteraction<T>(*iter,this->m_wall,m_k,iflag));
00078     this->m_inner_count+=(iflag ? 1 : 0);
00079   }
00080 
00081   console.XDebug() << "end CEWallInteractionGroup::Update()\n";
00082 }
00083 
00084 
00092 template<class T>
00093 void CEWallInteractionGroup<T>::applyForce(const Vec3& F)
00094 {
00095  
00096   int it=0;
00097   double d; // distance to move
00098   double df; // force difference
00099   double ef; // relative force error (df/|F|)
00100   Vec3 O_f=F.unit(); // direction of the applied force
00101   do{
00102     //std::cerr << "iteration: " << it << std::endl;
00103     // calculate local stiffness
00104     k_local=0.0;
00105     for(typename vector<CElasticWallInteraction<T> >::iterator iter=m_interactions.begin();
00106         iter!=m_interactions.end();
00107         iter++){      
00108       k_local+=iter->getStiffness();
00109     }
00110     //std::cerr << "local stiffness: " << k_local <<  std::endl;
00111     // get global K
00112     m_k_global=this->m_comm->sum_all(k_local);
00113     //std::cerr << "global stiffness: " << m_k_global <<  std::endl;
00114     if(m_k_global>0){
00115       // calculate local F
00116       Vec3 F_local=Vec3(0.0,0.0,0.0);
00117       for (
00118            typename vector<CElasticWallInteraction<T> >::iterator iter=m_interactions.begin();
00119            iter!=m_interactions.end();
00120            iter++
00121            ){
00122         if(iter->isInner()){
00123           Vec3 f_i=iter->getForce();
00124           F_local+=(f_i*O_f)*O_f; // add component of f_i in O_f direction
00125         }
00126       }
00127       //std::cerr << "local Force: " << F_local <<  std::endl;
00128       // get global F
00129       // by component (hack - fix later,i.e. sum_all for Vec3)
00130       double fgx=this->m_comm->sum_all(F_local.X());
00131       double fgy=this->m_comm->sum_all(F_local.Y());
00132       double fgz=this->m_comm->sum_all(F_local.Z());
00133       Vec3 F_global=Vec3(fgx,fgy,fgz);
00134       //std::cerr << "global Force: " << F_global <<  std::endl;
00135 
00136       // calc necessary wall movement
00137       df=(F+F_global)*O_f;
00138       d=df/m_k_global;
00139       ef=df/F.norm();
00140       //std::cerr << "d,df,ef: " << d  <<  " , " << df << " , " << ef << std::endl;
00141       // move the wall
00142       //std::cerr << " old wall pos: " <<  this->m_wall->getPos() <<  std::endl;
00143       this->m_wall->moveBy(d*O_f);
00144       //std::cerr << "new wall pos: " <<  this->m_wall->getPos() <<  std::endl;
00145       it++;
00146     } else {
00147       d=1e-5;
00148       ef=1;
00149       //std::cerr << "no contact" << std::endl;
00150       //std::cerr << "old wall pos: " <<  this->m_wall->getPos() <<  std::endl;
00151       this->m_wall->moveBy(d*O_f); 
00152       //std::cerr << "new wall pos: " <<  this->m_wall->getPos() <<  std::endl;
00153     } 
00154   } while((it<50)&&(ef>1e-3)); // check for convergence
00155 }
00156 
00157 
00158 template<class T>
00159 ostream& operator<<(ostream& ost,const CEWallInteractionGroup<T>& IG)
00160 {
00161   ost << "CEWallInteractionGroup" << endl << flush;
00162   ost << *(IG.m_wall) << endl << flush;
00163  
00164   return ost;
00165 }