ESyS-Particle  4.0.1
EWallInteraction.hpp
1 
2 // //
3 // Copyright (c) 2003-2011 by The University of Queensland //
4 // Earth Systems Science Computational Centre (ESSCC) //
5 // http://www.uq.edu.au/esscc //
6 // //
7 // Primary Business: Brisbane, Queensland, Australia //
8 // Licensed under the Open Software License version 3.0 //
9 // http://www.opensource.org/licenses/osl-3.0.php //
10 // //
12 
13 
14 #ifndef MODEL_EWALLINTERACTION_HPP
15 #define MODEL_EWALLINTERACTION_HPP
16 
25 template <class T>
27  AWallInteraction<T>(p,w,iflag)
28 {
29  double f=1.0;//this->m_p->getRad();;
30  // scale elastic param
31  if(!CParticle::getDo2dCalculations()){
32 // f*=this->m_p->getRad();
33  f*=3.141592654*this->m_p->getRad();
34  }
35  m_k=f*k;
36 }
37 
41 template <class T>
43 {
44 
45  double dist=(this->m_p->getPos()-this->m_wall->getOrigin())*this->m_wall->getNormal();
46  // console.XDebug() << "pos, rad, dist: " << this->m_p->getPos() << " " << this->m_p->getRad() << " " << dist << "\n";
47 
48  if(dist<this->m_p->getRad()){
49  Vec3 force=m_k*(this->m_p->getRad()-dist)*this->m_wall->getNormal();
50  Vec3 pos=this->m_p->getPos()-dist*this->m_wall->getNormal();
51 
52  /* Vec3 disp = this->m_p->getPos()- this->m_p->getOldPos() ;
53  Vec3 x_d = Vec3(1.0,0.0,0.0);
54  Vec3 y_d = Vec3(0.0,1.0,0.0);
55  Vec3 fy_vert = -m_k*disp.Y()*y_d;
56  Vec3 fx_vert = -m_k*disp.X()*x_d;
57  force = force+ fx_vert + fy_vert ;*/
58 
59 
60 // friction wall
61 /*
62  Vec3 x_d = Vec3(1.0,0.0,0.0);
63  Vec3 y_d = Vec3(0.0,1.0,0.0);
64  Vec3 fx_vert, fy_vert ;
65  double miu = 1.0;
66 
67  if ( this->m_p->getVel().X() >0.0) fx_vert = - miu*force.norm()*x_d;
68  else fx_vert = miu*force.norm()*x_d;
69 
70  if ( this->m_p->getVel().Y() >0.0) fy_vert = - miu*force.norm()*y_d;
71  else fy_vert = miu*force.norm()*y_d;
72 
73  force = force+ fx_vert + fy_vert ;
74  */
75  this->m_p->applyForce(force,pos);
76  if(this->m_inner_flag) this->m_wall->addForce(-1.0*force);
77  }
78 }
79 
83 template <class T>
85 {
86  Vec3 force=Vec3(0.0,0.0,0.0);
87  double dist=(this->m_p->getPos()-this->m_wall->getOrigin())*this->m_wall->getNormal();
88  if(dist<this->m_p->getRad()){
89  force=m_k*(this->m_p->getRad()-dist)*this->m_wall->getNormal();
90  }
91 
92  return -1.0*force;
93 }
94 
99 template <class T>
101 {
102  double res=0.0;
103  double dist=(this->m_p->getPos()-this->m_wall->getOrigin())*this->m_wall->getNormal();
104  if(dist<this->m_p->getRad()){
105  res=m_k;
106  }
107 
108  return res;
109 }
110 
111 
112 #endif