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 // 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 }