NDDEM
Contacts.h
Go to the documentation of this file.
1 
5 #ifndef CONTACTS
6 #define CONTACTS
7 
8 #include <cstdlib>
9 #include <cmath>
10 #include <cstdio>
11 #include <vector>
12 #include <cmath>
13 #include <map>
14 #include "Typedefs.h"
15 #include "Parameters.h"
16 #include "Tools.h"
17 #include "ContactList.h"
18 
21 template <int d>
22 class Contacts
23 {
24 public:
25  Contacts (Parameters<d> &PP) ;
26 
27  //int N ; ///< Number of particles // These and following should really be const but that's a pain for assignements
28  Parameters<d> *P ;
29  vector < double > Torquei ;
30  vector < double > Torquej ;
31  vector < double > vrel ;
32 
33  void particle_particle ( cv1d & Xi, cv1d & Vi, cv1d &Omegai, double ri, double mi,
34  cv1d & Xj, cv1d & Vj, cv1d &Omegaj, double rj, double mj, cp<d> & Contact, bool isdumptime) ;
35  void particle_wall ( cv1d & Vi, cv1d &Omegai, double ri, double mi,
36  cv1d & cn, cp<d> & Contact) ;
37  void particle_movingwall ( cv1d & Vi, cv1d & Omegai, double ri, double mi,
38  cv1d & cn, cv1d & Vj, cp<d> & Contact) ;
39  void particle_mesh ( cv1d & Xi, cv1d & Vi, cv1d &Omegai, double ri, double mi, cpm<d> & Contact) ;
40  void (Contacts::*particle_ghost) (cv1d & Xi, cv1d & Vi, cv1d &Omegai, double ri, double mi,
41  cv1d & Xj, cv1d & Vj, cv1d &Omegaj, double rj, double mj, cp<d> & Contact, bool isdumptime) ;
42 
43  void particle_ghost_regular (cv1d & Xi, cv1d & Vi, cv1d &Omegai, double ri, double mi,
44  cv1d & Xj, cv1d & Vj, cv1d &Omegaj, double rj, double mj, cp<d> & Contact, bool isdumptime)
45  {
46  vector <double> loc (d, 0) ;
47  loc=Xj ;
48  uint32_t gh=Contact.ghost, ghd=Contact.ghostdir ;
49  for (int n=0 ; gh>0 ; gh>>=1, ghd>>=1, n++)
50  {
51  if (gh&1)
52  loc[n] += P->Boundaries[n].delta * ((ghd&1)?-1:1) ;
53  }
54  return (particle_particle (Xi, Vi, Omegai, ri, mi, loc, Vj, Omegaj, rj, mj, Contact, isdumptime) ) ;
55  }
56 
57 
58 
59  void particle_ghost_LE (cv1d & Xi, cv1d & Vi, cv1d &Omegai, double ri, double mi,
60  cv1d & Xj, cv1d & Vj, cv1d &Omegaj, double rj, double mj, cp<d> & Contact, bool isdumptime) /*, FILE * debug = nullptr*/
61  {
62  vector <double> loc (d, 0) ;
63  //if (debug == nullptr) n= 1 ;
64  if ( P->Boundaries[0].Type != WallType::PBC_LE || (Contact.ghost & 1)==0)
65  {
66  loc=Xj ;
67  compute_normalpbcloc(loc, 0, Contact.ghost, Contact.ghostdir) ;
68  //printf(">> %g %g %g %g %g %g\n ", Xi[0], Xi[1], Xi[2], loc[0], loc[1], loc[2]) ;
69  particle_particle (Xi, Vi, Omegai, ri, mi, loc, Vj, Omegaj, rj, mj, Contact, isdumptime) ;
70  //if (debug != nullptr)
71  // fprintf(debug, "%g %g %g %g %g 0 %d\n", loc[0], loc[1], Vj[0], Vj[1], rj, Contact.j) ;
72  }
73  else
74  {
75  uint32_t gh=Contact.ghost, ghd=Contact.ghostdir ; // Handle pbc in first dim
76  vector <double> vel (d, 0) ; vel=Vj ;
77  loc=Xj ;
78  loc[0] += P->Boundaries[0].delta * ((ghd&1)?-1:1) ;
79  loc[1] += (ghd&1?-1:1)*P->Boundaries[0].displacement ;
80  double additionaldelta = 0 ;
81  if (loc[1] > P->Boundaries[1].xmax) {additionaldelta = -P->Boundaries[1].delta ;}
82  if (loc[1] < P->Boundaries[1].xmin) {additionaldelta = P->Boundaries[1].delta ;}
83  loc[1] += additionaldelta ;
84  vel[1] += (ghd&1?-1:1) * P->Boundaries[0].vel * P->Boundaries[0].delta ;
85 
86  gh>>=1 ; ghd>>=1 ;
87  compute_normalpbcloc (loc, 1, gh, ghd) ;
88  //printf(">> %g %g %g %g %g %g\n ", Xi[0], Xi[1], Xi[2], loc[0], loc[1], loc[2]) ;
89  particle_particle (Xi, Vi, Omegai, ri, mi, loc, vel, Omegaj, rj, mj, Contact, isdumptime) ;
90  }
91  }
92 
93  void particle_ghost_LE_cells (cv1d & Xi, cv1d & Vi, cv1d &Omegai, double ri, double mi,
94  cv1d & Xj, cv1d & Vj, cv1d &Omegaj, double rj, double mj, cp<d> & Contact, bool isdumptime) /*, FILE * debug = nullptr*/
95  {
96  vector <double> loc (d, 0) ;
97  //if (debug == nullptr) n= 1 ;
98 
99  if ( P->Boundaries[0].Type != WallType::PBC_LE || (Contact.ghost & 1)==0)
100  {
101  loc=Xj ;
102  compute_normalpbcloc(loc, 0, Contact.ghost, Contact.ghostdir) ;
103  particle_particle (Xi, Vi, Omegai, ri, mi, loc, Vj, Omegaj, rj, mj, Contact, isdumptime) ;
104  }
105  else
106  {
107  uint32_t gh=Contact.ghost, ghd=Contact.ghostdir ; // Handle pbc in first dim
108  vector <double> vel (d, 0) ; vel=Vj ;
109  loc=Xj ;
110  loc[0] += P->Boundaries[0].delta * ((ghd&1)?-1:1) ;
111  loc[1] += (ghd&1?-1:1)*P->Boundaries[0].displacement ;
112  vel[1] += (ghd&1?-1:1) * P->Boundaries[0].vel * P->Boundaries[0].delta ;
113 
114  gh>>=1 ; ghd>>=1 ;
115  compute_normalpbcloc (loc, 1, gh, ghd) ;
116  particle_particle (Xi, Vi, Omegai, ri, mi, loc, vel, Omegaj, rj, mj, Contact, isdumptime) ;
117  }
118  }
119 
120  void compute_normalpbcloc (v1d & loc, int startn, uint32_t gh, uint32_t ghd)
121  {
122  for (int n=startn ; gh>0 ; gh>>=1, ghd>>=1, n++)
123  {
124  if (gh&1)
125  loc[n] += P->Boundaries[n].delta * ((ghd&1)?-1:1) ;
126  }
127  }
128 
129  Action<d> Act ;
130 
131 private:
133  double contactlength, ovlp, dotvrelcn, Coulomb, tsprn, tsprcn ;
134  vector <double> cn, rri, rrj, vn, vt, Fn, Ft, tvec, Ftc, tspr, tsprc ;
135 
136 } ;
137 
138 
139 /*****************************************************************************************************
140  * *
141  * *
142  * *
143  * IMPLEMENTATIONS *
144  * *
145  * *
146  * *
147  * ***************************************************************************************************/
148 
149 template <int d>
150 Contacts<d>::Contacts (Parameters<d> &PP) : P(&PP)
151 {
152 
153  Torquei.resize(d*(d-1)/2, 0) ;
154  Torquej.resize(d*(d-1)/2, 0) ;
155  vrel.resize(d,0) ;
156  cn.resize(d,0) ;
157  rri.resize(d,0) ; rrj.resize(d,0) ; vn.resize(d,0) ; vt.resize(d,0) ;
158  Fn.resize(d,0) ; Ft.resize(d,0) ; tvec.resize(d,0) ; Ftc.resize(d,0) ;
159  tspr.resize(d,0) ; tsprc.resize(d,0) ;
160 
161  if (PP.contact_strategy == ContactStrategies::CELLS)
162  particle_ghost = &Contacts::particle_ghost_LE_cells ;
163  else
164  particle_ghost = &Contacts::particle_ghost_LE ;
165 
166 }
167 
168 //--------------------------------------------------------------------------------------
169 //---------------------- particle particle contact ----------------------------
170 template <int d>
171 void Contacts<d>::particle_particle (cv1d & Xi, cv1d & Vi, cv1d & Omegai, double ri, double mi,
172  cv1d & Xj, cv1d & Vj, cv1d & Omegaj, double rj, double mj, cp<d> & Contact, bool isdumptime)
173 {
174  double kn, kt, gamman, gammat ;
175  contactlength=Contact.contactlength ;
176 
177  ovlp=ri+rj-contactlength ;
178  if (ovlp<=0) {Act.setzero() ; Tools<d>::setzero(Contact.tspr) ; return ;}
179  //printf("%g %g %g %g\n", ri, rj, Xi[2], Xj[2]) ; fflush(stdout) ;
180  Tools<d>::vMinus(cn, Xi, Xj) ; //cn=(Xi-Xj) ;
181  cn /= contactlength ;
182 
183  //Relative velocity at contact
184  Tools<d>::vMul(rri, cn, ovlp/2.-ri) ; // rri = -cn * (ri-ovlp/2.) ;
185  Tools<d>::vMul(rrj, cn, rj-ovlp/2.) ; // rrj = cn * (rj-ovlp/2.) ;
186  vrel= Vi - Tools<d>::skewmatvecmult(Omegai, rri) - (Vj - Tools<d>::skewmatvecmult(Omegaj, rrj)) ; //TODO
187  Tools<d>::vMul (vn, cn, Tools<d>::dot(vrel,cn)) ; //vn=cn * (Tools<d>::dot(vrel,cn)) ;
188  Tools<d>::vMinus(vt, vrel, vn) ; //vt= vrel - vn ;
189 
190  Act.setvel(vn, vt) ;
191 
192  if (P->ContactModel == HERTZ)
193  {
194  double Rstar = (ri*rj)/(ri+rj) ;
195  double Mstar = (mi*mj)/(mi+mj) ;
196  double sqrtovlpR = sqrt(ovlp*Rstar) ;
197  kn = P->Kn * sqrtovlpR;
198  kt = P->Kt * sqrtovlpR;
199  gamman = P->Gamman * sqrt(Mstar * kn);
200  gammat = P->Gammat * sqrt(Mstar * kt);
201  }
202  else
203  {
204  kn=P->Kn ;
205  kt=P->Kt ;
206  gamman = P->Gamman ;
207  gammat = P->Gammat ;
208  }
209 
210  //Normal force
211  //Fn=cn*(ovlp*kn) - vn*gamman ;
212  Tools<d>::vMul(Fn, cn, ovlp*kn) ;
213  Tools<d>::vSubScaled (Fn, gamman, vn) ; //Fn=cn*(ovlp*kn) - vn*gamman ;
214  if (isdumptime)
215  {
216  Act.Fn_el = cn*(ovlp*kn) ;
217  Act.Fn_visc = - vn*gamman ;
218  }
219 
220 
221  //Tangential force computation: retrieve contact or create new contact
222  tspr=Contact.tspr ;
223  if (tspr.size()==0) tspr.resize(d,0) ;
224 
225  Tools<d>::vAddScaled (tspr, P->dt, vt) ; //tspr += vt*dt ;
226  Tools<d>::vSubScaled(tspr, Tools<d>::dot(tspr,cn), cn) ; // tspr -= cn * Tools<d>::dot(tspr,cn) ; //WARNING: might need an additional scaling so that |tsprnew|=|tspr|
227  Tools<d>::vMul(Ft, tspr, - kt) ; //Ft= tspr*(-Kt) ;
228  Coulomb=P->Mu*Tools<d>::norm(Fn) ;
229 
230  if (Tools<d>::norm(Ft) >= Coulomb)
231  {
232  if (Tools<d>::norm(tspr)>0)
233  {
234  Tools<d>::vMul(tvec, Ft, 1/Tools<d>::norm(Ft)) ; //tvec=Ft * (1/Tools<d>::norm(Ft)) ;
235  Tools<d>::vMul(Ftc, tvec, Coulomb) ; //Ftc = tvec * Coulomb ;
236  Ft=Ftc ;
237  Tools<d>::vMul(tspr, Ftc, -1/ kt) ; //tspr=Ftc*(-1/Kt) ;
238  }
239  else
240  Tools<d>::setzero(Ft) ;
241 
242  if (isdumptime) {Act.Ft_fric = Ft ; Act.Ft_isfric=true ; }
243  }
244  else
245  {
246  if (isdumptime) {Act.Ft_el = Ft ; Act.Ft_visc=-vt*gammat ; Act.Ft_isfric=false ; }
247  Tools<d>::vSubScaled(Ft, gammat, vt) ; //Ft -= (vt*Gammat) ;
248  }
249  Tools<d>::wedgeproduct(Torquei, rri, Ft) ;
250  Tools<d>::wedgeproduct(Torquej, Ft, rrj) ; //TODO check the minus sign Torquej = rrj^-Ft = Ft^rrj <= better for speed
251 
252  Contact.tspr=tspr ;
253  //Act.set_rev(-Fn, -Ft, Torquej) ;
254  //Act.set_fwd(Fn, Ft, Torquei) ;
255  Act.set(Fn, Ft, Torquei, Torquej) ;
256  return ;
257 }
258 
259 //---------------------- particle wall contact ----------------------------
260 template <int d>
261 void Contacts<d>::particle_wall ( cv1d & Vi, cv1d &Omegai, double ri, double mi,
262  cv1d & cn, cp<d> & Contact)
263  //int j, int orient, cp & Contact)
264 {
265  contactlength=Contact.contactlength ;
266  ovlp=ri-contactlength ;
267  if (ovlp<=0) {Act.setzero() ; Tools<d>::setzero(Contact.tspr) ; return ;}
268 
269  //Relative velocity at contact
270  rri = -cn * (ri-ovlp/2.) ;
271  vrel= Vi - Tools<d>::skewmatvecmult(Omegai, rri) ;
272  vn = cn * (Tools<d>::dot(vrel, cn)) ;
273  vt= vrel - vn ;
274  Act.setvel(vn, vt) ;
275 
276  double kn, kt, gamman, gammat ;
277  if (P->ContactModel == HERTZ)
278  {
279  double sqrtovlpR = sqrt(ovlp*ri) ;
280  kn = P->Kn * sqrtovlpR;
281  kt = P->Kt * sqrtovlpR;
282  gamman = P->Gamman * mi * sqrtovlpR;
283  gammat = P->Gammat * mi * sqrtovlpR;
284  }
285  else
286  {
287  kn=P->Kn ;
288  kt=P->Kt ;
289  gamman = P->Gamman ;
290  gammat = P->Gammat ;
291  }
292 
293  //Fn=cn*(ovlp*kn) - vn*gamman;
294  Tools<d>::vMul(Fn, cn, ovlp*kn) ;
295  Tools<d>::vSubScaled (Fn, gamman, vn) ; //Fn=cn*(ovlp*kn) - vn*gamman ;
296 
297  //Tangential force computation: retrieve contact or create new contact
298  //history=History[make_pair(i,-(2*j+k+1))] ;
299  tspr=Contact.tspr ; //history.second ;
300  if (tspr.size()==0) tspr.resize(d,0) ;
301 
302  tspr += vt*P->dt ;
303  tspr -= cn * Tools<d>::dot(tspr,cn) ; //WARNING: might need an additional scaling so that |tsprnew|=|tspr| Actually does not seem to change anything ...
304  Ft= tspr*(- kt) ;
305  Coulomb=P->Mu_wall*Tools<d>::norm(Fn) ;
306 
307  if (Tools<d>::norm(Ft) > Coulomb)
308  {
309  if (Tools<d>::norm(tspr)>0)
310  {
311  tvec=Ft * (1/Tools<d>::norm(Ft)) ;
312  Ftc = tvec * Coulomb ;
313  Ft=Ftc ;
314  tspr=Ftc*(-1/ kt) ;
315  }
316  else
317  Tools<d>::setzero(Ft) ;
318  }
319  else
320  Ft -= (vt*gammat) ;
321 
322  Torquei=Tools<d>::wedgeproduct(rri, Ft) ;
323 
324  Contact.tspr=tspr ;
325  Act.set(Fn, Ft, Torquei, Torquej) ;
326  return ;
327 }
328 
329 //---------------------- particle moving wall contact ----------------------------
330 template <int d>
331 void Contacts<d>::particle_movingwall ( cv1d & Vi, cv1d & Omegai, double ri, double mi,
332  cv1d & cn, cv1d & Vj, cp<d> & Contact)
333 {
334  contactlength=Contact.contactlength ;
335 
336  ovlp=ri-contactlength ;
337  if (ovlp<=0) {Act.setzero() ; return ;}
338 
339  //Relative velocity at contact
340  Tools<d>::vMul(rri, cn, ovlp/2.-ri) ; // rri = -cn * (ri-ovlp/2.) ;
341  vrel= Vi - Tools<d>::skewmatvecmult(Omegai, rri) - Vj ;
342  Tools<d>::vMul (vn, cn, Tools<d>::dot(vrel,cn)) ; //vn=cn * (Tools<d>::dot(vrel,cn)) ;
343  Tools<d>::vMinus(vt, vrel, vn) ; //vt= vrel - vn ;
344  Act.setvel(vn, vt) ;
345 
346  double kn, kt, gamman, gammat ;
347  if (P->ContactModel == HERTZ)
348  {
349  double sqrtovlpR = sqrt(ovlp*ri) ;
350  kn = P->Kn * sqrtovlpR;
351  kt = P->Kt * sqrtovlpR;
352  gamman = P->Gamman * mi * sqrtovlpR;
353  gammat = P->Gammat * mi * sqrtovlpR;
354  }
355  else
356  {
357  kn=P->Kn ;
358  kt=P->Kt ;
359  gamman = P->Gamman ;
360  gammat = P->Gammat ;
361  }
362 
363  //Normal force
364  //Fn=cn*(ovlp*kn) - vn*gamman ;
365  Tools<d>::vMul(Fn, cn, ovlp*kn) ;
366  Tools<d>::vSubScaled (Fn, gamman, vn) ; //Fn=cn*(ovlp*kn) - vn*gamman ;
367 
368  //Tangential force computation: retrieve contact or create new contact
369  tspr=Contact.tspr ;
370  if (tspr.size()==0) tspr.resize(d,0) ;
371 
372  Tools<d>::vAddScaled (tspr, P->dt, vt) ; //tspr += vt*dt ;
373  Tools<d>::vSubScaled(tspr, Tools<d>::dot(tspr,cn), cn) ; // tspr -= cn * Tools<d>::dot(tspr,cn) ; //WARNING: might need an additional scaling so that |tsprnew|=|tspr|
374  Tools<d>::vMul(Ft, tspr, - kt) ; //Ft= tspr*(-Kt) ;
375  Coulomb=P->Mu_wall*Tools<d>::norm(Fn) ;
376 
377  if (Tools<d>::norm(Ft) >= Coulomb)
378  {
379  if (Tools<d>::norm(tspr)>0)
380  {
381  Tools<d>::vMul(tvec, Ft, 1/Tools<d>::norm(Ft)) ; //tvec=Ft * (1/Tools<d>::norm(Ft)) ;
382  Tools<d>::vMul(Ftc, tvec, Coulomb) ; //Ftc = tvec * Coulomb ;
383  Ft=Ftc ;
384  Tools<d>::vMul(tspr, Ftc, -1/ kt) ; //tspr=Ftc*(-1/Kt) ;
385  }
386  else
387  Tools<d>::setzero(Ft) ;
388  }
389  else
390  Tools<d>::vSubScaled(Ft, gammat, vt) ; //Ft -= (vt*Gammat) ;
391 
392  Tools<d>::wedgeproduct(Torquei, rri, Ft) ;
393 
394  Contact.tspr=tspr ;
395  Act.set(Fn, Ft, Torquei, Torquej) ;
396  return ;
397 }
398 
399 //------------------------- Particle mesh contact --------------------------------------
400 template <int d>
401 void Contacts<d>::particle_mesh ( cv1d & Xi, cv1d & Vi, cv1d &Omegai, double ri, double mi, cpm<d> & Contact)
402 {
403  contactlength=Contact.contactlength ;
404 
405  ovlp=ri-contactlength ;
406  if (ovlp<=0) {Act.setzero() ; return ;}
407  Tools<d>::vMinus(cn, Xi, Contact.contactpoint) ; //cn=(Xi-Xj) ;
408  cn /= contactlength ;
409 
410  //Relative velocity at contact
411  Tools<d>::vMul(rri, cn, ovlp/2.-ri) ; // rri = -cn * (ri-ovlp/2.) ;
412  vrel= Vi - Tools<d>::skewmatvecmult(Omegai, rri) ;
413  Tools<d>::vMul (vn, cn, Tools<d>::dot(vrel,cn)) ; //vn=cn * (Tools<d>::dot(vrel,cn)) ;
414  Tools<d>::vMinus(vt, vrel, vn) ; //vt= vrel - vn ;
415  Act.setvel(vn, vt) ;
416 
417  double kn, kt, gamman, gammat ;
418  if (P->ContactModel == HERTZ)
419  {
420  double sqrtovlpR = sqrt(ovlp*ri) ;
421  kn = P->Kn * sqrtovlpR;
422  kt = P->Kt * sqrtovlpR;
423  gamman = P->Gamman * mi * sqrtovlpR;
424  gammat = P->Gammat * mi * sqrtovlpR;
425  }
426  else
427  {
428  kn=P->Kn ;
429  kt=P->Kt ;
430  gamman = P->Gamman ;
431  gammat = P->Gammat ;
432  }
433 
434  //Normal force
435  //Fn=cn*(ovlp*kn) - vn*gamman ; //TODO
436  Tools<d>::vMul(Fn, cn, ovlp*kn) ;
437  Tools<d>::vSubScaled (Fn, gamman, vn) ; //Fn=cn*(ovlp*kn) - vn*gamman ;
438 
439  //Tangential force computation: retrieve contact or create new contact
440  //history=History[make_pair(i,-(2*j+k+1))] ;
441  tspr=Contact.tspr ; //history.second ;
442  if (tspr.size()==0) tspr.resize(d,0) ;
443 
444  tspr += vt*P->dt ;
445  tspr -= cn * Tools<d>::dot(tspr,cn) ; //WARNING: might need an additional scaling so that |tsprnew|=|tspr| Actually does not seem to change anything ...
446  Ft= tspr*(- kt) ;
447  Coulomb=P->Mu_wall*Tools<d>::norm(Fn) ;
448 
449  if (Tools<d>::norm(Ft) > Coulomb)
450  {
451  if (Tools<d>::norm(tspr)>0)
452  {
453  tvec=Ft * (1/Tools<d>::norm(Ft)) ;
454  Ftc = tvec * Coulomb ;
455  Ft=Ftc ;
456  tspr=Ftc*(-1/ kt) ;
457  }
458  else
459  Tools<d>::setzero(Ft) ;
460  }
461  else
462  Ft -= (vt*gammat) ;
463 
464  Torquei=Tools<d>::wedgeproduct(rri, Ft) ;
465 
466  Contact.tspr=tspr ;
467  Act.set(Fn, Ft, Torquei, Torquej) ;
468  return ;
469 }
470 
471 
472 
473 #endif
Action
Handle force and torque contact information.
Definition: ContactList.h:19
Contacts
Calculate contact forces.
Definition: Contacts.h:23
Contacts::tsprc
vector< double > tsprc
Definition: Contacts.h:134
Contacts::particle_ghost_regular
void particle_ghost_regular(cv1d &Xi, cv1d &Vi, cv1d &Omegai, double ri, double mi, cv1d &Xj, cv1d &Vj, cv1d &Omegaj, double rj, double mj, cp< d > &Contact, bool isdumptime)
Calculate the particle to regular (non Lees-Edward) ghost contact.
Definition: Contacts.h:43
Contacts::vn
vector< double > vn
Definition: Contacts.h:134
Contacts::tspr
vector< double > tspr
Definition: Contacts.h:134
Contacts::vrel
vector< double > vrel
Relative velocity.
Definition: Contacts.h:31
Contacts::particle_ghost
void(Contacts::* particle_ghost)(cv1d &Xi, cv1d &Vi, cv1d &Omegai, double ri, double mi, cv1d &Xj, cv1d &Vj, cv1d &Omegaj, double rj, double mj, cp< d > &Contact, bool isdumptime)
Function pointer to the function to calculate the ghost-to-particle contact forces.
Definition: Contacts.h:40
Contacts::Ft
vector< double > Ft
Definition: Contacts.h:134
Contacts::vt
vector< double > vt
Definition: Contacts.h:134
Contacts::Ftc
vector< double > Ftc
Definition: Contacts.h:134
Contacts::rri
vector< double > rri
Definition: Contacts.h:134
Contacts::P
Parameters< d > * P
Pointer to the Parameters structure.
Definition: Contacts.h:28
Contacts::ovlp
double ovlp
Definition: Contacts.h:133
Contacts::Coulomb
double Coulomb
Definition: Contacts.h:133
Contacts::tsprn
double tsprn
Definition: Contacts.h:133
Contacts::Torquej
vector< double > Torquej
Torque on particle j.
Definition: Contacts.h:30
Contacts::dotvrelcn
double dotvrelcn
Definition: Contacts.h:133
Contacts::contactlength
double contactlength
Definition: Contacts.h:133
Contacts::compute_normalpbcloc
void compute_normalpbcloc(v1d &loc, int startn, uint32_t gh, uint32_t ghd)
Move ghosts through the regular periodic boundary conditions (non Lees-Edward).
Definition: Contacts.h:120
Contacts::tvec
vector< double > tvec
Definition: Contacts.h:134
Contacts::cn
vector< double > cn
Definition: Contacts.h:134
Contacts::Fn
vector< double > Fn
Definition: Contacts.h:134
Contacts::particle_ghost_LE
void particle_ghost_LE(cv1d &Xi, cv1d &Vi, cv1d &Omegai, double ri, double mi, cv1d &Xj, cv1d &Vj, cv1d &Omegaj, double rj, double mj, cp< d > &Contact, bool isdumptime)
Calculate the particle to (potentially) Lees-Edward ghost contact.
Definition: Contacts.h:59
Contacts::tsprcn
double tsprcn
Definition: Contacts.h:133
Contacts::Act
Action< d > Act
Resulting Action.
Definition: Contacts.h:129
Contacts::Torquei
vector< double > Torquei
Torque on particle i.
Definition: Contacts.h:29
Contacts::rrj
vector< double > rrj
Definition: Contacts.h:134
Contacts::particle_ghost_LE_cells
void particle_ghost_LE_cells(cv1d &Xi, cv1d &Vi, cv1d &Omegai, double ri, double mi, cv1d &Xj, cv1d &Vj, cv1d &Omegaj, double rj, double mj, cp< d > &Contact, bool isdumptime)
Calculate the particle to (potentially) Lees-Edward ghost contact for the cell contact finding (sligh...
Definition: Contacts.h:93
Parameters
Generic class to handle the simulation set up.
Definition: Parameters.h:128
Parameters::contact_strategy
ContactStrategies contact_strategy
Strategy for the contact detection.
Definition: Parameters.h:188
Tools
Static class to handle multi-dimensional mathematics, and more. It gets specialised for speed with te...
Definition: Tools.h:101
Tools::vAddScaled
static void vAddScaled(v1d &res, double v, cv1d &a, cv1d &b)
Addition of 3 vectors in-place.
Definition: Tools.h:183
Tools::vMul
static void vMul(v1d &res, cv1d &a, double b)
Multiply a vector by a scalar in-place.
Definition: Tools.h:175
Tools::vSubScaled
static void vSubScaled(v1d &res, double v, cv1d &a)
Subtraction of a scaled vector .
Definition: Tools.h:186
Tools::vMinus
static void vMinus(v1d &res, cv1d &a, cv1d &b)
Difference of 2 vectors in-place .
Definition: Tools.h:177
cp
Contact properties class.
Definition: ContactList.h:69
cpm
Contact properties for mesh contacts (mainly), including contact point location, specialising cp.
Definition: ContactList.h:177
Contact
@ Contact
Definition: Coarsing.h:59
Tools::setzero
static void setzero(v2d &a)
Set a matrix to zero in-place.
Definition: Tools.h:117
Tools::wedgeproduct
static v1d wedgeproduct(cv1d &a, cv1d &b)
Wedge product of vectors.
Definition: Tools.h:710
Tools::dot
static double dot(cv1d &a, cv1d &b)
Dot product.
Definition: Tools.h:127
Contacts::particle_particle
void particle_particle(cv1d &Xi, cv1d &Vi, cv1d &Omegai, double ri, double mi, cv1d &Xj, cv1d &Vj, cv1d &Omegaj, double rj, double mj, cp< d > &Contact, bool isdumptime)
Force & torque between 2 particles.
Definition: Contacts.h:171
vel
@ vel
Definition: Typedefs.h:19
Contacts::particle_mesh
void particle_mesh(cv1d &Xi, cv1d &Vi, cv1d &Omegai, double ri, double mi, cpm< d > &Contact)
Force & torque between particle and mesh.
Definition: Contacts.h:401
Contacts::particle_movingwall
void particle_movingwall(cv1d &Vi, cv1d &Omegai, double ri, double mi, cv1d &cn, cv1d &Vj, cp< d > &Contact)
Force & torque between a particle and a moving wall. Vj is the velocity of the wall at the contact po...
Definition: Contacts.h:331
Contacts::Contacts
Contacts(Parameters< d > &PP)
Definition: Contacts.h:150
Tools::skewmatvecmult
static v1d skewmatvecmult(cv1d &M, cv1d &v)
Multiply the skew symetric matrix M with vector v.
Definition: Tools.h:569
cv1d
const vector< double > cv1d
Definition: Typedefs.h:13
Contacts::particle_wall
void particle_wall(cv1d &Vi, cv1d &Omegai, double ri, double mi, cv1d &cn, cp< d > &Contact)
Force & torque between a particle and a wall.
Definition: Contacts.h:261
v1d
vector< double > v1d
Definition: Typedefs.h:9
Tools::norm
static double norm(const vector< double > &a)
Norm of a vector.
Definition: Tools.h:119
CELLS
@ CELLS
Definition: Typedefs.h:22
HERTZ
@ HERTZ
Definition: Typedefs.h:21
WallType::PBC_LE
@ PBC_LE
d
uint d
uint32_t
unsigned int uint32_t
Definition: stdint.h:126
Generated by doxygen 1.9.1