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 #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