ESyS-Particle  4.0.1
Damping.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_DAMPING_HPP
00015 #define MODEL_DAMPING_HPP
00016 
00017 #include "Model/Damping.h"
00018 //#include "Foundation/console.h"
00019 
00020 template <class T>
00021 double CDamping<T>::s_limit2=1e-12; // default error limit : 1e-6 
00022 template <class T>
00023 int CDamping<T>::s_flops = 0;
00024 
00025 
00035 template <class T>
00036 CDamping<T>::CDamping(T* P,const Vec3& V,double visc,double dt,int mi)
00037 {
00038   m_p=P;
00039   m_vref=V;
00040   m_visc=visc;
00041   m_dt=dt;
00042   m_maxiter=mi;
00043   m_force=Vec3(0.0,0.0,0.0);
00044 }
00045 
00052 template <class T>
00053 CDamping<T>::CDamping(T* P,const CDampingIGP& param)
00054 {
00055   m_p=P;
00056   m_vref=param.getVRef();
00057   m_visc=param.getVisc();
00058   m_dt=param.getTimeStep();
00059   m_maxiter=param.getMaxIter();
00060   m_force=Vec3(0.0,0.0,0.0);
00061 }
00062 
00069 template <class T>
00070 CDamping<T>::CDamping(T* P,CDampingIGP* param)
00071 {
00072   m_p=P;
00073   m_vref=param->getVRef();
00074   m_visc=param->getVisc();
00075   m_dt=param->getTimeStep();
00076   m_maxiter=param->getMaxIter();
00077   m_force=Vec3(0.0,0.0,0.0);
00078 }
00079 
00083 template <class T>
00084 CDamping<T>::~CDamping()
00085 {}
00086 
00087 template <class T>
00088 void CDamping<T>::setTimeStepSize(double dt)
00089 {
00090   m_dt = dt;
00091 }
00092 
00098 template <class T>
00099 void CDamping<T>::calcForces()
00100 {
00101   m_E_diss=0.0;// zero dissipated energy
00102   Vec3 v=m_p->getVel();
00103   Vec3 v_rel=Vec3(0.0,0.0,0.0);
00104   Vec3 frc=m_p->getForce();
00105   
00106   double s=m_p->getInvMass();     
00107   double mass=m_p->getMass();
00108   double error=1.0;
00109   int count=0;
00110   while((error*error>s_limit2) & (count<m_maxiter)){ // 1 flop
00111     count++;
00112     Vec3 v_old=v_rel;
00113     v_rel=v-m_vref+s*m_dt*(frc-v_rel*m_visc*mass);        // 16 flops
00114     error=(v_rel-v_old).norm2();                     // 8 flops  
00115   }
00116   if(count>=m_maxiter){
00117     //console.Warning() << "damping diverges for particle  " << m_p->getID() << "error= " << error << "\n";
00118     v_rel=Vec3(0.0,0.0,0.0);
00119   }
00120   m_force=-1.0*m_visc*v_rel*mass;
00121   m_p->applyForce(m_force,m_p->getPos());  // 3 flops
00122   // m_E_diss=0.5*(v*v-(v_rel+m_vref)*(v_rel+m_vref))*m_p->getMass();
00123   m_E_diss=m_visc*v_rel*v*m_dt;
00124 }
00125 
00131 template <class T>
00132 typename CDamping<T>::ScalarFieldFunction CDamping<T>::getScalarFieldFunction(const string& name)
00133 {
00134   typename CDamping<T>::ScalarFieldFunction sf;
00135 
00136   if(name=="dissipated_energy"){
00137     sf=&CDamping<T>::getDissipatedEnergy;
00138   } else {
00139     sf=NULL;
00140     cerr << "ERROR - invalid name for interaction scalar  access function" << endl;
00141   }
00142 
00143   return sf;
00144 }
00145 
00151 template <class T>
00152 typename CDamping<T>::CheckedScalarFieldFunction CDamping<T>::getCheckedScalarFieldFunction(const string& name)
00153 {
00154   typename CDamping<T>::CheckedScalarFieldFunction sf;
00155 
00156   sf=NULL;
00157   cerr << "ERROR - invalid name for interaction scalar  access function" << endl;
00158   
00159   return sf;
00160 }
00161 
00162 
00168 template <class T>
00169 typename CDamping<T>::VectorFieldFunction CDamping<T>::getVectorFieldFunction(const string& name)
00170 {
00171   typename CDamping<T>::VectorFieldFunction vf;
00172 
00173   if (name=="force"){
00174     vf=&CDamping<T>::getForce;
00175   } else {
00176       vf=NULL;
00177       cerr << "ERROR - invalid name for interaction vector access function" << endl;
00178   }
00179   
00180   return vf;
00181 }
00182 
00186 template <class T>
00187 double CDamping<T>::getDissipatedEnergy() const
00188 {
00189   return m_E_diss;
00190 }
00191 
00192 template <class T>
00193 Vec3 CDamping<T>::getForce() const
00194 {
00195   return m_force;
00196 }
00197 
00204 template <class T>
00205 bool CDamping<T>::hasTag(int tag,int mask) const
00206 {
00207  int tag1=m_p->getTag();
00208 
00209   return ((tag1 & mask)==(tag & mask));
00210 }
00211 
00215 template <class T>
00216 vector<int> CDamping<T>::getAllID() const
00217 {
00218   vector<int> res;
00219 
00220   res.push_back(m_p->getID());
00221 
00222   return res;
00223 }
00224 
00225 #endif