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