ESyS-Particle  4.0.1
RotDamping.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 #ifndef __ROT_DAMPING_HPP
00014 #define __ROT_DAMPING_HPP
00015 
00016 template <class T>
00017 double CRotDamping<T>::s_limit2=1e-12; // default error limit : 1e-6 
00018 template <class T>
00019 int CRotDamping<T>::s_flops = 0;
00020 
00021 
00028 template <class T>
00029 CRotDamping<T>::CRotDamping(T* P,CDampingIGP* param)
00030 {
00031   m_p=P;
00032   m_vref=param->getVRef();
00033   m_visc=param->getVisc();
00034   m_dt=param->getTimeStep();
00035   m_maxiter=param->getMaxIter();
00036 }
00037 
00041 template <class T>
00042 CRotDamping<T>::~CRotDamping()
00043 {
00044 }
00045 
00046 template <class T>
00047 void CRotDamping<T>::setTimeStepSize(double dt)
00048 {
00049   m_dt = dt;
00050 }
00051 
00057 template <class T>
00058 void CRotDamping<T>::calcForces()
00059 {
00060   m_E_diss=0.0;// zero dissipated energy
00061   Vec3 v=m_p->getAngVel();
00062   Vec3 v_rel=Vec3(0.0,0.0,0.0);
00063   Vec3 frc=m_p->getMoment();
00064   
00065   double s=1.0/m_p->getInertRot();       
00066   double in=m_p->getInertRot();
00067   double mass=m_p->getMass();
00068 
00069   double error=1.0;
00070   int count=0;
00071   while((error*error>s_limit2) & (count<m_maxiter)){ // 1 flop
00072     count++;
00073     Vec3 v_old=v_rel;
00074     v_rel=v-m_vref+s*m_dt*(frc-v_rel*m_visc*in);        // 16 flops
00075     error=(v_rel-v_old).norm2();                     // 8 flops  
00076   }
00077   if(count>=m_maxiter){
00078     //console.Warning() << "damping diverges for particle  " << m_p->getID() << "error= " << error << "\n";
00079     v_rel=Vec3(0.0,0.0,0.0);
00080   }
00081   m_force=-1.0*m_visc*v_rel*mass;
00082   m_p->applyMoment(-1.0*m_visc*v_rel*in);  // 3 flops
00083   m_E_diss=m_visc*v_rel*v*m_dt; // wrong
00084 }
00085 
00091 template <class T>
00092 typename CRotDamping<T>::ScalarFieldFunction CRotDamping<T>::getScalarFieldFunction(const string& name)
00093 {
00094   typename CRotDamping<T>::ScalarFieldFunction sf;
00095 
00096   if(name=="dissipated_energy"){
00097     sf=&CRotDamping<T>::getDissipatedEnergy;
00098   } else {
00099     sf=NULL;
00100     cerr << "ERROR - invalid name for interaction scalar  access function" << endl;
00101   }
00102 
00103   return sf;
00104 }
00105 
00111 template <class T>
00112 typename CRotDamping<T>::CheckedScalarFieldFunction CRotDamping<T>::getCheckedScalarFieldFunction(const string& name)
00113 {
00114   typename CRotDamping<T>::CheckedScalarFieldFunction sf;
00115 
00116   sf=NULL;
00117   cerr << "ERROR - invalid name for interaction scalar  access function" << endl;
00118   
00119   return sf;
00120 }
00121 
00122 
00128 template <class T>
00129 typename CRotDamping<T>::VectorFieldFunction CRotDamping<T>::getVectorFieldFunction(const string& name)
00130 {
00131   typename CRotDamping<T>::VectorFieldFunction vf;
00132 
00133   if (name=="force"){
00134     vf=&CRotDamping<T>::getForce;
00135   } else {
00136     vf=NULL;
00137     cerr << "ERROR - invalid name for interaction vector access function" << endl;
00138   }
00139   
00140   return vf;
00141 }
00142 
00146 template <class T>
00147 double CRotDamping<T>::getDissipatedEnergy() const
00148 {
00149   return m_E_diss;
00150 }
00151 
00152 template <class T>
00153 Vec3 CRotDamping<T>::getForce() const
00154 {
00155   return m_force;
00156 }
00157 
00164 template <class T>
00165 bool CRotDamping<T>::hasTag(int tag,int mask) const
00166 {
00167  int tag1=m_p->getTag();
00168 
00169   return ((tag1 & mask)==(tag & mask));
00170 }
00171 
00175 template <class T>
00176 vector<int> CRotDamping<T>::getAllID() const
00177 {
00178   vector<int> res;
00179 
00180   res.push_back(m_p->getID());
00181 
00182   return res;
00183 }
00184 
00185 #endif // __ROT_DAMPING_HPP