ESyS-Particle  4.0.1
RotDamping.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 #ifndef __ROT_DAMPING_HPP
14 #define __ROT_DAMPING_HPP
15 
16 template <class T>
17 double CRotDamping<T>::s_limit2=1e-12; // default error limit : 1e-6
18 template <class T>
20 
21 
28 template <class T>
30 {
31  m_p=P;
32  m_vref=param->getVRef();
33  m_visc=param->getVisc();
34  m_dt=param->getTimeStep();
35  m_maxiter=param->getMaxIter();
36 }
37 
41 template <class T>
43 {
44 }
45 
46 template <class T>
47 void CRotDamping<T>::setTimeStepSize(double dt)
48 {
49  m_dt = dt;
50 }
51 
57 template <class T>
59 {
60  m_E_diss=0.0;// zero dissipated energy
61  Vec3 v=m_p->getAngVel();
62  Vec3 v_rel=Vec3(0.0,0.0,0.0);
63  Vec3 frc=m_p->getMoment();
64 
65  double s=1.0/m_p->getInertRot();
66  double in=m_p->getInertRot();
67  double mass=m_p->getMass();
68 
69  double error=1.0;
70  int count=0;
71  while((error*error>s_limit2) & (count<m_maxiter)){ // 1 flop
72  count++;
73  Vec3 v_old=v_rel;
74  v_rel=v-m_vref+s*m_dt*(frc-v_rel*m_visc*in); // 16 flops
75  error=(v_rel-v_old).norm2(); // 8 flops
76  }
77  if(count>=m_maxiter){
78  //console.Warning() << "damping diverges for particle " << m_p->getID() << "error= " << error << "\n";
79  v_rel=Vec3(0.0,0.0,0.0);
80  }
81  m_force=-1.0*m_visc*v_rel*mass;
82  m_p->applyMoment(-1.0*m_visc*v_rel*in); // 3 flops
83  m_E_diss=m_visc*v_rel*v*m_dt; // wrong
84 }
85 
91 template <class T>
93 {
95 
96  if(name=="dissipated_energy"){
98  } else {
99  sf=NULL;
100  cerr << "ERROR - invalid name for interaction scalar access function" << endl;
101  }
102 
103  return sf;
104 }
105 
111 template <class T>
113 {
115 
116  sf=NULL;
117  cerr << "ERROR - invalid name for interaction scalar access function" << endl;
118 
119  return sf;
120 }
121 
122 
128 template <class T>
130 {
132 
133  if (name=="force"){
135  } else {
136  vf=NULL;
137  cerr << "ERROR - invalid name for interaction vector access function" << endl;
138  }
139 
140  return vf;
141 }
142 
146 template <class T>
148 {
149  return m_E_diss;
150 }
151 
152 template <class T>
154 {
155  return m_force;
156 }
157 
164 template <class T>
165 bool CRotDamping<T>::hasTag(int tag,int mask) const
166 {
167  int tag1=m_p->getTag();
168 
169  return ((tag1 & mask)==(tag & mask));
170 }
171 
175 template <class T>
176 vector<int> CRotDamping<T>::getAllID() const
177 {
178  vector<int> res;
179 
180  res.push_back(m_p->getID());
181 
182  return res;
183 }
184 
185 #endif // __ROT_DAMPING_HPP