ESyS-Particle  4.0.1
brokenEWallInteractionGroup.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   double 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     if(iflag){
00080       if(!CParticle::getDo2dCalculations()){
00081         k_local+=m_k*((*iter)->getRad()); // update local K
00082       } else {
00083         k_local+=m_k;
00084       }
00085     }
00086   }
00087   // get global K
00088   m_k_global=this->m_comm->sum_all(k_local);
00089 
00090   console.XDebug() << "end CEWallInteractionGroup::Update()\n";
00091 }
00092 
00093 
00101 template<class T>
00102 void CEWallInteractionGroup<T>::applyForce(const Vec3& F)
00103 {
00104   int it=0;
00105   double d;
00106   Vec3 O_f=F.unit(); // direction of the applied force
00107   do{
00108     // calculate local F
00109     Vec3 F_local=Vec3(0.0,0.0,0.0);
00110     for (
00111       typename vector<CElasticWallInteraction<T> >::iterator iter=m_interactions.begin();
00112       iter!=m_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)/m_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 
00135 
00136 template<class T>
00137 ostream& operator<<(ostream& ost,const CEWallInteractionGroup<T>& IG)
00138 {
00139   ost << "CEWallInteractionGroup" << endl << flush;
00140   ost << *(IG.m_wall) << endl << flush;
00141  
00142   return ost;
00143 }