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