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