ESyS-Particle  4.0.1
EWallInteraction.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 #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