ESyS-Particle  2.3.2
Damping.hpp
Go to the documentation of this file.
1 
2 // //
3 // Copyright (c) 2003-2014 by The University of Queensland //
4 // Centre for Geoscience Computing //
5 // http://earth.uq.edu.au/centre-geoscience-computing //
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 
14 #ifndef MODEL_DAMPING_HPP
15 #define MODEL_DAMPING_HPP
16 
17 #include "Model/Damping.h"
18 //#include "Foundation/console.h"
19 
20 template <class T>
21 double CDamping<T>::s_limit2=1e-12; // default error limit : 1e-6
22 template <class T>
23 int CDamping<T>::s_flops = 0;
24 
25 
35 template <class T>
36 CDamping<T>::CDamping(T* P,const Vec3& V,double visc,double dt,int mi)
37 {
38  m_p=P;
39  m_vref=V;
40  m_visc=visc;
41  m_dt=dt;
42  m_maxiter=mi;
43  m_force=Vec3(0.0,0.0,0.0);
44 }
45 
52 template <class T>
54 {
55  m_p=P;
56  m_vref=param.getVRef();
57  m_visc=param.getVisc();
58  m_dt=param.getTimeStep();
59  m_maxiter=param.getMaxIter();
60  m_force=Vec3(0.0,0.0,0.0);
61 }
62 
69 template <class T>
71 {
72  m_p=P;
73  m_vref=param->getVRef();
74  m_visc=param->getVisc();
75  m_dt=param->getTimeStep();
76  m_maxiter=param->getMaxIter();
77  m_force=Vec3(0.0,0.0,0.0);
78 }
79 
83 template <class T>
85 {}
86 
87 template <class T>
89 {
90  m_dt = dt;
91 }
92 
98 template <class T>
100 {
101  m_E_diss=0.0;// zero dissipated energy
102  Vec3 v=m_p->getVel();
103  Vec3 v_rel=Vec3(0.0,0.0,0.0);
104  Vec3 frc=m_p->getForce();
105 
106  double s=m_p->getInvMass();
107  double mass=m_p->getMass();
108  double error=1.0;
109  int count=0;
110  while((error*error>s_limit2) & (count<m_maxiter)){ // 1 flop
111  count++;
112  Vec3 v_old=v_rel;
113  v_rel=v-m_vref+s*m_dt*(frc-v_rel*m_visc*mass); // 16 flops
114  error=(v_rel-v_old).norm2(); // 8 flops
115  }
116  if(count>=m_maxiter){
117  //console.Warning() << "damping diverges for particle " << m_p->getID() << "error= " << error << "\n";
118  v_rel=Vec3(0.0,0.0,0.0);
119  }
120  m_force=-1.0*m_visc*v_rel*mass;
121  m_p->applyForce(m_force,m_p->getPos()); // 3 flops
122  // m_E_diss=0.5*(v*v-(v_rel+m_vref)*(v_rel+m_vref))*m_p->getMass();
123  m_E_diss=m_visc*v_rel*v*m_dt;
124 }
125 
131 template <class T>
133 {
135 
136  if(name=="dissipated_energy"){
138  } else {
139  sf=NULL;
140  cerr << "ERROR - invalid name for interaction scalar access function" << endl;
141  }
142 
143  return sf;
144 }
145 
151 template <class T>
153 {
155 
156  sf=NULL;
157  cerr << "ERROR - invalid name for interaction scalar access function" << endl;
158 
159  return sf;
160 }
161 
162 
168 template <class T>
170 {
172 
173  if (name=="force"){
175  } else {
176  vf=NULL;
177  cerr << "ERROR - invalid name for interaction vector access function" << endl;
178  }
179 
180  return vf;
181 }
182 
186 template <class T>
188 {
189  return m_E_diss;
190 }
191 
192 template <class T>
194 {
195  return m_force;
196 }
197 
204 template <class T>
205 bool CDamping<T>::hasTag(int tag,int mask) const
206 {
207  int tag1=m_p->getTag();
208 
209  return ((tag1 & mask)==(tag & mask));
210 }
211 
215 template <class T>
216 vector<int> CDamping<T>::getAllID() const
217 {
218  vector<int> res;
219 
220  res.push_back(m_p->getID());
221 
222  return res;
223 }
224 
225 #endif