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 #ifndef MODEL_EWALLINTERACTION_HPP 00015 #define MODEL_EWALLINTERACTION_HPP 00016 00025 template <class T> 00026 CElasticWallInteraction<T>::CElasticWallInteraction(T* p,CWall* w,double k,bool iflag): 00027 AWallInteraction<T>(p,w,iflag) 00028 { 00029 double f=1.0;//this->m_p->getRad();; 00030 // scale elastic param 00031 if(!CParticle::getDo2dCalculations()){ 00032 // f*=this->m_p->getRad(); 00033 f*=3.141592654*this->m_p->getRad(); 00034 } 00035 m_k=f*k; 00036 } 00037 00041 template <class T> 00042 void CElasticWallInteraction<T>::calcForces() 00043 { 00044 00045 double dist=(this->m_p->getPos()-this->m_wall->getOrigin())*this->m_wall->getNormal(); 00046 // console.XDebug() << "pos, rad, dist: " << this->m_p->getPos() << " " << this->m_p->getRad() << " " << dist << "\n"; 00047 00048 if(dist<this->m_p->getRad()){ 00049 Vec3 force=m_k*(this->m_p->getRad()-dist)*this->m_wall->getNormal(); 00050 Vec3 pos=this->m_p->getPos()-dist*this->m_wall->getNormal(); 00051 00052 /* Vec3 disp = this->m_p->getPos()- this->m_p->getOldPos() ; 00053 Vec3 x_d = Vec3(1.0,0.0,0.0); 00054 Vec3 y_d = Vec3(0.0,1.0,0.0); 00055 Vec3 fy_vert = -m_k*disp.Y()*y_d; 00056 Vec3 fx_vert = -m_k*disp.X()*x_d; 00057 force = force+ fx_vert + fy_vert ;*/ 00058 00059 00060 // friction wall 00061 /* 00062 Vec3 x_d = Vec3(1.0,0.0,0.0); 00063 Vec3 y_d = Vec3(0.0,1.0,0.0); 00064 Vec3 fx_vert, fy_vert ; 00065 double miu = 1.0; 00066 00067 if ( this->m_p->getVel().X() >0.0) fx_vert = - miu*force.norm()*x_d; 00068 else fx_vert = miu*force.norm()*x_d; 00069 00070 if ( this->m_p->getVel().Y() >0.0) fy_vert = - miu*force.norm()*y_d; 00071 else fy_vert = miu*force.norm()*y_d; 00072 00073 force = force+ fx_vert + fy_vert ; 00074 */ 00075 this->m_p->applyForce(force,pos); 00076 if(this->m_inner_flag) this->m_wall->addForce(-1.0*force); 00077 } 00078 } 00079 00083 template <class T> 00084 Vec3 CElasticWallInteraction<T>::getForce() 00085 { 00086 Vec3 force=Vec3(0.0,0.0,0.0); 00087 double dist=(this->m_p->getPos()-this->m_wall->getOrigin())*this->m_wall->getNormal(); 00088 if(dist<this->m_p->getRad()){ 00089 force=m_k*(this->m_p->getRad()-dist)*this->m_wall->getNormal(); 00090 } 00091 00092 return -1.0*force; 00093 } 00094 00099 template <class T> 00100 double CElasticWallInteraction<T>::getStiffness() 00101 { 00102 double res=0.0; 00103 double dist=(this->m_p->getPos()-this->m_wall->getOrigin())*this->m_wall->getNormal(); 00104 if(dist<this->m_p->getRad()){ 00105 res=m_k; 00106 } 00107 00108 return res; 00109 } 00110 00111 00112 #endif