MercuryDPM  0.10
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StatisticsPoint.hcc
Go to the documentation of this file.
1 //Copyright (c) 2013-2014, The MercuryDPM Developers Team. All rights reserved.
2 //For the list of developers, see <http://www.MercuryDPM.org/Team>.
3 //
4 //Redistribution and use in source and binary forms, with or without
5 //modification, are permitted provided that the following conditions are met:
6 // * Redistributions of source code must retain the above copyright
7 // notice, this list of conditions and the following disclaimer.
8 // * Redistributions in binary form must reproduce the above copyright
9 // notice, this list of conditions and the following disclaimer in the
10 // documentation and/or other materials provided with the distribution.
11 // * Neither the name MercuryDPM nor the
12 // names of its contributors may be used to endorse or promote products
13 // derived from this software without specific prior written permission.
14 //
15 //THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
16 //ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
17 //WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
18 //DISCLAIMED. IN NO EVENT SHALL THE MERCURYDPM DEVELOPERS TEAM BE LIABLE FOR ANY
19 //DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
20 //(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
21 //LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
22 //ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
23 //(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
24 //SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
25 
26 std::ostream& operator<<(std::ostream& os, const CG S) {
27  if (S==HeavisideSphere) os << "HeavisideSphere";
28  else if (S==Gaussian) os << "Gaussian";
29  else if (S==Polynomial) os << "Polynomial";
30  return os;
31 }
32 
33 template <StatType T>
34 StatisticsVector<T>* StatisticsPoint<T>::gb;
35 
36 template <StatType T>
37 void StatisticsPoint<T>::set_zero() {
38  Nu = 0.0;
39  Density = 0.0;
40  Momentum.set_zero();
41  DisplacementMomentum.set_zero();
42  Displacement.set_zero();
43  MomentumFlux.set_zero();
44  DisplacementMomentumFlux.set_zero();
45  EnergyFlux.set_zero();
46  TangentialStress.set_zero();
47  NormalStress.set_zero();
48  TangentialStress.set_zero();
49  NormalTraction.set_zero();
50  TangentialTraction.set_zero();
51  Fabric.set_zero();
52  CollisionalHeatFlux.set_zero();
53  Dissipation = 0.0;
54  Potential = 0.0;
55  LocalAngularMomentum.set_zero();
56  LocalAngularMomentumFlux.set_zero();
57  ContactCoupleStress.set_zero();
58 }
59 
60 template <StatType T>
61 StatisticsPoint<T> StatisticsPoint<T>::getSquared() {
62  StatisticsPoint<T> P;
63  P.Nu = sqr(Nu);
64  P.Density = sqr(Density);
65  P.Momentum = square(Momentum);
66  P.DisplacementMomentum = square(DisplacementMomentum);
67  P.Displacement = square(Displacement);
68  P.MomentumFlux = square(MomentumFlux);
69  P.DisplacementMomentumFlux = square(DisplacementMomentumFlux);
70  P.EnergyFlux = square(EnergyFlux);
71  P.NormalStress = square(NormalStress);
72  P.TangentialStress = square(TangentialStress);
73  P.NormalTraction = square(NormalTraction);
74  P.TangentialTraction = square(TangentialTraction);
75  P.Fabric = square(Fabric);
76  P.CollisionalHeatFlux = square(CollisionalHeatFlux);
77  P.Dissipation = sqr(Dissipation);
78  P.Potential = sqr(Potential);
79  LocalAngularMomentum = square(P.LocalAngularMomentum);
80  LocalAngularMomentumFlux = square(P.LocalAngularMomentumFlux);
81  ContactCoupleStress = square(P.ContactCoupleStress);
82  return P;
83 }
84 
85 template <StatType T>
86 StatisticsPoint<T>& StatisticsPoint<T>::operator=(const StatisticsPoint<T> &P)
87 {
88  Nu = P.Nu;
89  Density = P.Density;
90  Momentum = P.Momentum;
91  DisplacementMomentum = P.DisplacementMomentum;
92  Displacement = P.Displacement;
93  MomentumFlux = P.MomentumFlux;
94  DisplacementMomentumFlux = P.DisplacementMomentumFlux;
95  EnergyFlux = P.EnergyFlux;
96  NormalStress = P.NormalStress;
97  TangentialStress = P.TangentialStress;
98  NormalTraction = P.NormalTraction;
99  TangentialTraction = P.TangentialTraction;
100  Fabric = P.Fabric;
101  CollisionalHeatFlux = P.CollisionalHeatFlux;
102  Dissipation = P.Dissipation;
103  Potential = P.Potential;
104  LocalAngularMomentum = P.LocalAngularMomentum;
105  LocalAngularMomentumFlux = P.LocalAngularMomentumFlux;
106  ContactCoupleStress = P.ContactCoupleStress;
107  return *this;
108 }
109 
110 template <StatType T>
111 StatisticsPoint<T>& StatisticsPoint<T>::operator+=(const StatisticsPoint<T> &P)
112 {
113  Nu += P.Nu;
114  Density += P.Density;
115  Momentum += P.Momentum;
116  DisplacementMomentum += P.DisplacementMomentum;
117  Displacement += P.Displacement;
118  MomentumFlux += P.MomentumFlux;
119  DisplacementMomentumFlux += P.DisplacementMomentumFlux;
120  EnergyFlux += P.EnergyFlux;
121  NormalStress += P.NormalStress;
122  TangentialStress += P.TangentialStress;
123  NormalTraction += P.NormalTraction;
124  TangentialTraction += P.TangentialTraction;
125  Fabric += P.Fabric;
126  CollisionalHeatFlux += P.CollisionalHeatFlux;
127  Dissipation += P.Dissipation;
128  Potential += P.Potential;
129  LocalAngularMomentum += P.LocalAngularMomentum;
130  LocalAngularMomentumFlux += P.LocalAngularMomentumFlux;
131  ContactCoupleStress += P.ContactCoupleStress;
132  return *this;
133 }
134 
135 template <StatType T>
136 StatisticsPoint<T>& StatisticsPoint<T>::operator-=(const StatisticsPoint<T> &P)
137 {
138  Nu -= P.Nu;
139  Density -= P.Density;
140  Momentum -= P.Momentum;
141  DisplacementMomentum -= P.DisplacementMomentum;
142  Displacement -= P.Displacement;
143  MomentumFlux -= P.MomentumFlux;
144  DisplacementMomentumFlux -= P.DisplacementMomentumFlux;
145  EnergyFlux -= P.EnergyFlux;
146  NormalStress -= P.NormalStress;
147  TangentialStress -= P.TangentialStress;
148  NormalTraction -= P.NormalTraction;
149  TangentialTraction -= P.TangentialTraction;
150  Fabric -= P.Fabric;
151  CollisionalHeatFlux -= P.CollisionalHeatFlux;
152  Dissipation -= P.Dissipation;
153  Potential -= P.Potential;
154  LocalAngularMomentum -= P.LocalAngularMomentum;
155  LocalAngularMomentumFlux -= P.LocalAngularMomentumFlux;
156  ContactCoupleStress -= P.ContactCoupleStress;
157  return *this;
158 }
159 
160 template <StatType T>
161 StatisticsPoint<T>& StatisticsPoint<T>::operator/=(const Mdouble a)
162 {
163  Nu /= a;
164  Density /= a;
165  Momentum /= a;
166  DisplacementMomentum /= a;
167  Displacement /= a;
168  MomentumFlux /= a;
169  DisplacementMomentumFlux /= a;
170  EnergyFlux /= a;
171  NormalStress /= a;
172  TangentialStress /= a;
173  NormalTraction /= a;
174  TangentialTraction /= a;
175  Fabric /= a;
176  CollisionalHeatFlux /= a;
177  Dissipation /= a;
178  Potential /= a;
179  LocalAngularMomentum /= a;
180  LocalAngularMomentumFlux /= a;
181  ContactCoupleStress /= a;
182  return *this;
183 }
184 
185 //in the first average, one timestep in Displacement is missing
186 template <StatType T>
187 void StatisticsPoint<T>::firstTimeAverage(const int n) {
188  Nu /= n;
189  Density /= n;
190  Momentum /= n;
191  if (n==1) {
192  DisplacementMomentum.set_zero();
193  Displacement.set_zero();
194  DisplacementMomentumFlux.set_zero();
195  } else {
196  DisplacementMomentum /= n-1;
197  Displacement /= n-1;
198  DisplacementMomentumFlux /= n-1;
199  }
200  MomentumFlux /= n;
201  EnergyFlux /= n;
202  NormalStress /= n;
203  TangentialStress /= n;
204  NormalTraction /= n;
205  TangentialTraction /= n;
206  Fabric /= n;
207  CollisionalHeatFlux /= n;
208  Dissipation /= n;
209  Potential /= n;
210  LocalAngularMomentum /= n;
211  LocalAngularMomentumFlux /= n;
212  ContactCoupleStress /= n;
213 }
214 
215 template <StatType T>
216 Mdouble StatisticsPoint<T>::get_distance2(const Vec3D &P) {
217  return sqr(P.X-Position.X)+sqr(P.Y-Position.Y)+sqr(P.Z-Position.Z);
218 }
219 
220 template<StatType T>
221 Vec3D StatisticsPoint<T>::cross(Vec3D P, Vec3D &Q) {
222  return Vec3D(P.Y*Q.Z-P.Z*Q.Y, P.Z*Q.X-P.X*Q.Z, P.X*Q.Y-P.Y*Q.X);
223 }
224 
225 template<StatType T> Matrix3D StatisticsPoint<T>::MatrixCross(Vec3D P, Matrix3D &Q) {
226  return Matrix3D(
227  P.Y*Q.ZX-P.Z*Q.YX, P.Z*Q.XX-P.X*Q.ZX, P.X*Q.YX-P.Y*Q.XX,
228  P.Y*Q.ZY-P.Z*Q.YY, P.Z*Q.XY-P.X*Q.ZY, P.X*Q.YY-P.Y*Q.XY,
229  P.Y*Q.ZZ-P.Z*Q.YZ, P.Z*Q.XZ-P.X*Q.ZZ, P.X*Q.YZ-P.Y*Q.XZ);
230 }
231 
232 template <StatType T>
233 Mdouble StatisticsPoint<T>::dot(const Vec3D &P, const Vec3D &Q) {
234  return P.X*Q.X + P.Y*Q.Y + P.Z*Q.Z;
235 }
236 
237 template <StatType T>
238 Mdouble StatisticsPoint<T>::CG_function(const Vec3D &PI) {
239  Mdouble dist2 = get_distance2(PI);
240  if (get_CG_type()==HeavisideSphere) {
241  return (get_w2()<dist2)?0.0:get_CG_invvolume();
242  } else if (get_CG_type()==Gaussian) {
243  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume() * exp(-dist2/(2.0*get_w2()));
244  } else if (get_CG_type()==Polynomial) {
245  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume()*evaluatePolynomial(sqrt(dist2)/get_cutoff());
246  } else { std::cerr << "error in CG_function" << std::endl; exit(-1); }
247 }
248 
249 
250 template <StatType T>
251 Mdouble StatisticsPoint<T>::CG_function_2D(const Vec3D &PI) {
252  Mdouble dist2 = get_distance2(PI);
253  if (get_CG_type()==HeavisideSphere) {
254  if (get_w2()<dist2) return 0.0;
255  else {
256  Mdouble wn = sqrt(get_w2()-dist2);
257  return get_CG_invvolume()*2.0*wn;
259  //return get_CG_invvolume()*(std::min(wn,gb->getObjectDistanceLeft()) +std::min(wn,gb->getObjectDistanceRight()));
260  }
261  } else if (get_CG_type()==Gaussian) {
262  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume() * exp(-dist2/(2.0*get_w2()));
263  } else if (get_CG_type()==Polynomial) {
264  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume()*evaluatePolynomial(sqrt(dist2)/get_cutoff());
265  } else { std::cerr << "error in CG_function_2D" << std::endl; exit(-1); }
266 }
267 
268 template <StatType T>
269 Mdouble StatisticsPoint<T>::CG_function_1D(const Vec3D &PI) {
270  Mdouble dist2 = get_distance2(PI);
271  if (get_CG_type()==HeavisideSphere) {
272  return (get_w2()<dist2)?0.0:(get_CG_invvolume()*constants::pi*(get_w2()-dist2));
273  } else if (get_CG_type()==Gaussian) {
274  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume() * exp(-dist2/(2.0*get_w2()));
275  } else if (get_CG_type()==Polynomial) {
276  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume()*evaluatePolynomial(sqrt(dist2)/get_cutoff());
277  } else { std::cerr << "error in CG_function_1D" << std::endl; exit(-1); }
278 }
279 
280 template <StatType T>
281 Vec3D StatisticsPoint<T>::CG_gradient(const Vec3D &P, const Mdouble phi) {
282  //CG_type_todo
283  if (get_CG_type()==Gaussian) {
284  return (P-Position)*(phi/get_w2());
285  } else if (get_CG_type()==Polynomial) {
286  Mdouble r=GetLength(Position-P)/get_w();
287  return get_CG_invvolume()*evaluatePolynomialGradient(r)/r/get_w2()*(Position-P);
288  } else { std::cerr << "error in CG_gradient" << std::endl; exit(-1); }
289 }
290 
292 template <StatType T>
293 Vec3D StatisticsPoint<T>::CG_integral_gradient(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance) {
294  Vec3D P_P1 = Position - P1;
295  Vec3D P_P2 = Position - P2;
296  Mdouble a = dot(P_P1,P1_P2_normal);
297  Mdouble b = dot(P_P2,P1_P2_normal);
298  Vec3D tangential = P_P1 - a * P1_P2_normal;
299  //CG_type_todo
300  if (get_CG_type()==Gaussian) {
301  Mdouble wsq2 = sqrt(2*get_w2());
302  Mdouble f = sqrt(2*constants::pi*get_w2());
303  return Vec3D(0.0,0.0,(exp(-sqr(a/wsq2))-exp(-sqr(b/wsq2)))/f/(P2.Z-P1.Z));
304  } else if (get_CG_type()==Polynomial) {
305  double wn2 = get_w2() - dot(tangential,tangential);
306  if ((wn2<=0) | (a*fabs(a)>=wn2) | (b*fabs(b)<=-wn2)) {
307  return Vec3D(0,0,0);
308  } else {
310  double wn = sqrt(wn2);
311  double w = get_w();
312  double delta = w*1e-3;
313  double I = get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,tangential.GetLength()/w)*w / P1_P2_distance;
314  Vec3D delta_P_P1 = P_P1 +Vec3D(delta,0,0);
315  Vec3D delta_P_P2 = P_P2 +Vec3D(delta,0,0);
316  a = dot(delta_P_P1,P1_P2_normal);
317  b = dot(delta_P_P2,P1_P2_normal);
318  tangential = delta_P_P1 - a * P1_P2_normal;
319  wn = sqrt(get_w2() - dot(tangential,tangential));
320  double Ix = get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,tangential.GetLength()/w)*w / P1_P2_distance;
321  delta_P_P1 = P_P1 +Vec3D(0,delta,0);
322  delta_P_P2 = P_P2 +Vec3D(0,delta,0);
323  a = dot(delta_P_P1,P1_P2_normal);
324  b = dot(delta_P_P2,P1_P2_normal);
325  tangential = delta_P_P1 - a * P1_P2_normal;
326  wn = sqrt(get_w2() - dot(tangential,tangential));
327  double Iy = get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,tangential.GetLength()/w)*w / P1_P2_distance;
328  delta_P_P1 = P_P1 +Vec3D(0,0,delta);
329  delta_P_P2 = P_P2 +Vec3D(0,0,delta);
330  a = dot(delta_P_P1,P1_P2_normal);
331  b = dot(delta_P_P2,P1_P2_normal);
332  tangential = delta_P_P1 - a * P1_P2_normal;
333  wn = sqrt(get_w2() - dot(tangential,tangential));
334  double Iz = get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,tangential.GetLength()/w)*w / P1_P2_distance;
335  return Vec3D(Ix-I,Iy-I,Iz-I)/delta;
336  }
337  } else { return Vec3D(0,0,0);}
338  //~ if (get_CG_type()==Gaussian) return Vec3D(0.0,0.0,0.0);
339  //~ else { std::cerr << "error in CG_function" << std::endl; exit(-1); }
340 }
341 
342 template <StatType T>
343 Mdouble StatisticsPoint<T>::CG_integral_gradient_1D(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance) {
344  Vec3D P_P1 = Position - P1;
345  Vec3D P_P2 = Position - P2;
346  Mdouble a = dot(P_P1,P1_P2_normal);
347  Mdouble b = dot(P_P2,P1_P2_normal);
348  Vec3D tangential = P_P1 - a * P1_P2_normal;
349  //CG_type_todo
350  if (get_CG_type()==Gaussian) {
351  Mdouble wsq2 = sqrt(2*get_w2());
352  Mdouble f = sqrt(2*constants::pi*get_w2());
353  return (exp(-sqr(a/wsq2))-exp(-sqr(b/wsq2)))/f/(P2.Z-P1.Z);
354  } else if (get_CG_type()==Polynomial) {
355  double wn2 = get_w2() - dot(tangential,tangential);
356  if ((wn2<=0) | (a*fabs(a)>=wn2) | (b*fabs(b)<=-wn2)) {
357  return 0;
358  } else if ((P1_P2_distance<1e-12)|(GetLength2(tangential)<1e-24)) {
359  //if the normal is parallel to the averaging direction
360  std::cout << "normal is parallel to the averaging direction: "
361  << " P1_P2_distance " << P1_P2_distance
362  << " tangential " << tangential << std::endl;
364  return 0;//CG_gradient_1D(P1,CG_function(P1));
365  } else {
366  double wn = sqrt(wn2);
367  double w = get_w();
369  double delta = w*3e-6;
370  //double I = get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,dot(tangential,tangential)/w)*w / P1_P2_distance;
371  Vec3D delta_P_P1 = P_P1 +Vec3D(delta,delta,delta);
372  Vec3D delta_P_P2 = P_P2 +Vec3D(delta,delta,delta);
373  a = dot(delta_P_P1,P1_P2_normal);
374  b = dot(delta_P_P2,P1_P2_normal);
375  tangential = delta_P_P1 - a * P1_P2_normal;
376  wn = sqrt(get_w2() - dot(tangential,tangential));
377  double I2 = get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,dot(tangential,tangential)/w)*w / P1_P2_distance;
378 
379  delta_P_P1 = P_P1 -Vec3D(delta,delta,delta);
380  delta_P_P2 = P_P2 -Vec3D(delta,delta,delta);
381  a = dot(delta_P_P1,P1_P2_normal);
382  b = dot(delta_P_P2,P1_P2_normal);
383  tangential = delta_P_P1 - a * P1_P2_normal;
384  wn = sqrt(get_w2() - dot(tangential,tangential));
385  double I1 = get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,dot(tangential,tangential)/w)*w / P1_P2_distance;
386 
387  return (I2-I1)/(2.*delta);
388  }
389  } else return 0;
390 }
391 
398 template <StatType T>
400 Mdouble StatisticsPoint<T>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi UNUSED)
401 {
402  //apply the cutoff in normal direction:
403  Vec3D P_P1 = Position - P1;
404  Mdouble a = dot(P_P1,P1_P2_normal);
405  if ( a > get_cutoff() ) return 0.0;
406  Vec3D P_P2 = Position - P2;
407  Mdouble b = dot(P_P2,P1_P2_normal);
408  if ( -b > get_cutoff() ) return 0.0;
409  //apply the cutoff in tangential direction:
410  Vec3D tangential = P_P1 - a * P1_P2_normal;
411  if (dot(tangential,tangential)>get_cutoff2()) return 0.0;
412  //evaluate:
413  if (get_CG_type()==HeavisideSphere) {
414  Mdouble wn2 = get_w2() - dot(tangential,tangential);
415  if ((wn2<=0) | (a*fabs(a)>=wn2) | (b*fabs(b)<=-wn2)) {
416  return 0;
417  } else {
418  Mdouble wn = sqrt(wn2);
419  return get_CG_invvolume()*( std::min(b,wn)-std::max(a,-wn) ) / P1_P2_distance;
420  }
421  } else if (get_CG_type()==Gaussian) { //Gaussian
422  static Mdouble InvVolumeExp = compute_Gaussian_invvolume(gb->get_dim_particle()-1);
423  static Mdouble w_sqrt_2 = constants::sqrt_2*get_w();
424  static Mdouble P1_P2_distance_for_cutoff=-1, InvVolumeErf=-1;
425  if (P1_P2_distance_for_cutoff!=P1_P2_distance) {
426  P1_P2_distance_for_cutoff = P1_P2_distance;
427  Mdouble amax = -get_cutoff();
428  Mdouble bmax = get_cutoff()+P1_P2_distance;
429  InvVolumeErf = P1_P2_distance/(
430  erf(bmax/w_sqrt_2)*bmax+w_sqrt_2/constants::sqrt_pi*exp(-bmax*bmax/(w_sqrt_2*w_sqrt_2))
431  -erf(amax/w_sqrt_2)*amax-w_sqrt_2/constants::sqrt_pi*exp(-amax*amax/(w_sqrt_2*w_sqrt_2)));
432  // std::cout << "InvVolumeErf" << InvVolumeErf << " InvVolumeExp" << InvVolumeExp << " f" << get_CG_invvolume() * sqrt(2*constants::pi*get_w2()) << std::endl;
433  }
434  Mdouble psi=exp(-dot(tangential,tangential)/(2.0*get_w2())) * InvVolumeExp
435  * ( erf(b/w_sqrt_2) - erf(a/w_sqrt_2) ) * InvVolumeErf /2./P1_P2_distance;
436  rpsi=Position*psi;
437  return psi;
438  } else if (get_CG_type()==Polynomial) {
439  double wn2 = get_w2() - dot(tangential,tangential);
440  if ((wn2<=0) | (a*fabs(a)>=wn2) | (b*fabs(b)<=-wn2)) {
441  return 0;
442  } else {
443  double wn = sqrt(wn2);
444  double w = get_w();
445  return get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,tangential.GetLength()/w)*w / P1_P2_distance;
446  }
447  //CG_type_todo
448  } else { std::cerr << "error in CG_integral" << std::endl; exit(-1); }
449 }
450 
452 template <StatType T>
453 Mdouble StatisticsPoint<T>::CG_integral_2D(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Mdouble& rpsi_scalar)
454 {
455  //apply the cutoff in normal direction:
456  Vec3D P_P1 = Position - P1;
457  Mdouble a = dot(P_P1,P1_P2_normal);
458  if ( a > get_cutoff() ) return 0.0;
459  Vec3D P_P2 = Position - P2;
460  Mdouble b = dot(P_P2,P1_P2_normal);
461  if ( -b > get_cutoff() ) return 0.0;
462  //apply the cutoff in tangential direction:
463  Vec3D tangential = P_P1 - a * P1_P2_normal;
464  if (dot(tangential,tangential)>get_cutoff2()) return 0.0;
465  if (get_CG_type()==HeavisideSphere) {
466  Mdouble wn2 = get_w2() - dot(tangential,tangential);
467  if ((wn2<=0) || (a*fabs(a)>=wn2) || (b*fabs(b)<=-wn2)) {
468  return 0;
469  } else {
470  Mdouble wn = sqrt(wn2);
471  //if the normal is parallel to the averaging direction
472  if (std::max(fabs(a),fabs(b))<1e-20) {
473  return get_CG_invvolume() * 2 * wn;
474  }
475  return get_CG_invvolume() / P1_P2_distance
476  *( ((b>= wn)?(constants::pi*wn2/2.0):( b*sqrt(wn2-sqr(b))+wn2*asin(b/wn)))
477  +((a<=-wn)?(constants::pi*wn2/2.0):(-a*sqrt(wn2-sqr(a))-wn2*asin(a/wn))));
478  }
479  } else if (get_CG_type()==Gaussian) { //Gaussian
480  if (dot(P1_P2_normal,P1_P2_normal)<1e-20) {
481  //since we average parallel to the force line, the CG integral is shaped like the CG function
483  return get_CG_invvolume() * exp(-dot(P_P1,P_P1)/(2.0*get_w2()));
484  }
485  static Mdouble InvVolumeExp = compute_Gaussian_invvolume(gb->get_dim_particle()-2);
486  static Mdouble w_sqrt_2 = constants::sqrt_2*get_w();
487  static Mdouble P1_P2_distance_for_cutoff=-1, InvVolumeErf=-1;
488  if (P1_P2_distance_for_cutoff!=P1_P2_distance) {
489  P1_P2_distance_for_cutoff = P1_P2_distance;
490  Mdouble amax = -get_cutoff();
491  Mdouble bmax = get_cutoff()+P1_P2_distance;
492  InvVolumeErf = P1_P2_distance/(
493  erf(bmax/w_sqrt_2)*bmax+w_sqrt_2/constants::sqrt_pi*exp(-bmax*bmax/(w_sqrt_2*w_sqrt_2))
494  -erf(amax/w_sqrt_2)*amax-w_sqrt_2/constants::sqrt_pi*exp(-amax*amax/(w_sqrt_2*w_sqrt_2))
495  ) / averagingVolume();
496  //std::cout << "InvVolumeErf" << InvVolumeErf << " InvVolumeExp" << InvVolumeExp << " f" << get_CG_invvolume() * sqrt(2*constants::pi*get_w2()) << std::endl;
497  }
498  //Note: use dot(t,t), GetLength2(t), as dot only works on non-averaged directions
499  //we also define rpsi
500  Mdouble psi = exp(-dot(tangential,tangential)/(2.0*get_w2())) * InvVolumeExp
501  * ( erf(b/w_sqrt_2) - erf(a/w_sqrt_2) ) * InvVolumeErf /2./P1_P2_distance;
502  Mdouble phi1 = get_CG_invvolume() * exp(-dot(P_P1,P_P1)/(2.0*get_w2()));
503  Mdouble phi2 = get_CG_invvolume() * exp(-dot(P_P2,P_P2)/(2.0*get_w2()));
504  rpsi_scalar = -a/P1_P2_distance*psi + get_w2()/P1_P2_distance/P1_P2_distance*(phi1-phi2);
505  return psi;
506  } else if (get_CG_type()==Polynomial) {
507  double wn2 = get_w2() - dot(tangential,tangential);
508  if ((wn2<=0) | (a*fabs(a)>=wn2) | (b*fabs(b)<=-wn2)) {
509  return 0;
510  } else if (P1_P2_distance<1e-20) {
511  //if the normal is parallel to the averaging direction
512  return CG_function_1D(P_P1);
513  } else {
514  double wn = sqrt(wn2);
515  double w = get_w();
516  return get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,tangential.GetLength()/w)*w / P1_P2_distance;
517  }
518  //CG_type_todo
519  } else { std::cerr << "error in CG_integral_2D" << std::endl; exit(-1); }
520 }
521 
523 template <StatType T>
524 Mdouble StatisticsPoint<T>::CG_integral_1D(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Mdouble& rpsi_scalar)
525 {
526  //apply the cutoff in normal direction:
527  Vec3D P_P1 = Position - P1;
528  Mdouble a = dot(P_P1,P1_P2_normal);
529  if ( a > get_cutoff() ) return 0.0;
530  Vec3D P_P2 = Position - P2;
531  Mdouble b = dot(P_P2,P1_P2_normal);
532  if ( -b > get_cutoff() ) return 0.0;
533  //apply the cutoff in tangential direction:
534  Vec3D tangential = P_P1 - a * P1_P2_normal;
535  if (dot(tangential,tangential)>get_cutoff2()) return 0.0;
536  if (get_CG_type()==HeavisideSphere) {
537  Mdouble wn2 = get_w2() - dot(tangential,tangential);
538  if ((wn2<=0) | (a*fabs(a)>=wn2) | (b*fabs(b)<=-wn2)) {
539  return 0;
540  } else {
541  //if the normal is parallel to the averaging direction
542  if (std::max(fabs(a),fabs(b))<1e-20) {
543  return get_CG_invvolume() * constants::pi * wn2;
544  }
545  Mdouble wn = sqrt(wn2);
546  return get_CG_invvolume() / P1_P2_distance
547  *( ((b>= wn)?( 2.0/3.0*constants::pi*wn2*wn):(b*constants::pi*(wn2-sqr(b)/3.0)))
548  -((a<=-wn)?(-2.0/3.0*constants::pi*wn2*wn):(a*constants::pi*(wn2-sqr(a)/3.0))) );
549  }
550  } else if (get_CG_type()==Gaussian) { //Gaussian
551  if (fabs(b-a)<1e-20) {
552  return get_CG_invvolume() * exp(-dot(P_P1,P_P1)/(2.0*get_w2()));
553  }
554  static Mdouble w_sqrt_2 = constants::sqrt_2*get_w();
555  static Mdouble P1_P2_distance_for_cutoff=-1, InvVolumeErf=-1;
556  if (P1_P2_distance_for_cutoff!=P1_P2_distance) {
557  P1_P2_distance_for_cutoff = P1_P2_distance;
558  Mdouble amax = -get_cutoff();
559  Mdouble bmax = get_cutoff()+P1_P2_distance;
560  InvVolumeErf = P1_P2_distance/(
561  erf(bmax/w_sqrt_2)*bmax+w_sqrt_2/constants::sqrt_pi*exp(-bmax*bmax/(w_sqrt_2*w_sqrt_2))
562  -erf(amax/w_sqrt_2)*amax-w_sqrt_2/constants::sqrt_pi*exp(-amax*amax/(w_sqrt_2*w_sqrt_2))
563  ) / averagingVolume();
564  //std::cout << "InvVolumeErf" << InvVolumeErf << " InvVolumeExp" << InvVolumeExp << " f" << get_CG_invvolume() * sqrt(2*constants::pi*get_w2()) << std::endl;
565  }
566  //Note: use dot(t,t), GetLength2(t), as dot only works on non-averaged directions
567  Mdouble psi = ( erf(b/w_sqrt_2) - erf(a/w_sqrt_2) ) * InvVolumeErf /2./P1_P2_distance;
568  Mdouble phi1 = get_CG_invvolume() * exp(-dot(P_P1,P_P1)/(2.0*get_w2()));
569  Mdouble phi2 = get_CG_invvolume() * exp(-dot(P_P2,P_P2)/(2.0*get_w2()));
570  rpsi_scalar = -a/P1_P2_distance*psi + get_w2()/P1_P2_distance/P1_P2_distance*(phi1-phi2);
571  return psi;
572  } else if (get_CG_type()==Polynomial) {
573  //~ std::cout << setprecision(12) << P1_P2_normal << std::endl;
574  double wn2 = get_w2() - dot(tangential,tangential);
575  if ((wn2<=0) | (a*fabs(a)>=wn2) | (b*fabs(b)<=-wn2)) {
576  return 0;
577  } else if ((P1_P2_distance<1e-12)|(GetLength2(tangential)<1e-24)) {
578  //~ } else if (fabs(b-a)<1e-12) {
579  //if the normal is parallel to the averaging direction
580  //~ std::cout << "normal is parallel to the averaging direction: "
581  //~ << " P_P1 " << P1_P2_distance
582  //~ << " psi= " << CG_function_1D(P1) << std::endl;
583  return CG_function_1D(P1);
584  } else {
585  double wn = sqrt(wn2);
586  double w = get_w();
587  return get_CG_invvolume()*evaluateIntegral(std::max(a,-wn)/w,std::min(b,wn)/w,dot(tangential,tangential)/w)*w / P1_P2_distance;
588  }
589  } else { std::cerr << "error in CG_integral_1D" << std::endl; exit(-1); }
590  std::cout<<"eind testje"<< rpsi_scalar<<std::endl;
591 }
592 
593 template <StatType T>
594 std::string StatisticsPoint<T>::write_variable_names() {
595  std::stringstream ss;
596  ss<< "Nu "
597  << "Density "
598  << "MomentumX "
599  << "MomentumY "
600  << "MomentumZ "
601  << "DisplacementMomentumX "
602  << "DisplacementMomentumY "
603  << "DisplacementMomentumZ "
604  << "DisplacementXX "
605  << "DisplacementXY "
606  << "DisplacementXZ "
607  << "DisplacementYY "
608  << "DisplacementYZ "
609  << "DisplacementZZ "
610  << "MomentumFluxXX "
611  << "MomentumFluxXY "
612  << "MomentumFluxXZ "
613  << "MomentumFluxYY "
614  << "MomentumFluxYZ "
615  << "MomentumFluxZZ "
616  << "DisplacementMomentumFluxXX "
617  << "DisplacementMomentumFluxXY "
618  << "DisplacementMomentumFluxXZ "
619  << "DisplacementMomentumFluxYY "
620  << "DisplacementMomentumFluxYZ "
621  << "DisplacementMomentumFluxZZ "
622  << "EnergyFluxX "
623  << "EnergyFluxY "
624  << "EnergyFluxZ "
625  << "NormalStressXX "
626  << "NormalStressXY "
627  << "NormalStressXZ "
628  << "NormalStressYX "
629  << "NormalStressYY "
630  << "NormalStressYZ "
631  << "NormalStressZX "
632  << "NormalStressZY "
633  << "NormalStressZZ "
634  << "TangentialStressXX "
635  << "TangentialStressXY "
636  << "TangentialStressXZ "
637  << "TangentialStressYX "
638  << "TangentialStressYY "
639  << "TangentialStressYZ "
640  << "TangentialStressZX "
641  << "TangentialStressZY "
642  << "TangentialStressZZ "
643  << "NormalTractionX "
644  << "NormalTractionY "
645  << "NormalTractionZ "
646  << "TangentialTractionX "
647  << "TangentialTractionY "
648  << "TangentialTractionZ "
649  << "FabricXX "
650  << "FabricXY "
651  << "FabricXZ "
652  << "FabricYY "
653  << "FabricYZ "
654  << "FabricZZ "
655  << "CollisionalHeatFluxX "
656  << "CollisionalHeatFluxY "
657  << "CollisionalHeatFluxZ "
658  << "Dissipation "
659  << "Potential "
660  << "LocalAngularMomentumX "
661  << "LocalAngularMomentumY "
662  << "LocalAngularMomentumZ "
663  << "LocalAngularMomentumFluxXX "
664  << "LocalAngularMomentumFluxXY "
665  << "LocalAngularMomentumFluxXZ "
666  << "LocalAngularMomentumFluxYX "
667  << "LocalAngularMomentumFluxYY "
668  << "LocalAngularMomentumFluxYZ "
669  << "LocalAngularMomentumFluxZX "
670  << "LocalAngularMomentumFluxZY "
671  << "LocalAngularMomentumFluxZZ "
672  << "ContactCoupleStressXX "
673  << "ContactCoupleStressXY "
674  << "ContactCoupleStressXZ "
675  << "ContactCoupleStressYX "
676  << "ContactCoupleStressYY "
677  << "ContactCoupleStressYZ "
678  << "ContactCoupleStressZX "
679  << "ContactCoupleStressZY "
680  << "ContactCoupleStressZZ ";
681  return ss.str();
682 }
683 
684 template <StatType T>
685 std::string StatisticsPoint<T>::print() const
686 {
687  std::stringstream ss;
688  ss<< "Nu " << Nu
689  <<", Density " << Density
690  << ",\n Momentum " << Momentum
691  << ",\n DisplacementMomentum " << DisplacementMomentum
692  << ",\n Displacement " << Displacement
693  << ",\n MomentumFlux " << MomentumFlux
694  << ",\n DisplacementMomentumFlux " << DisplacementMomentumFlux
695  << ",\n EnergyFlux " << EnergyFlux
696  << ",\n NormalStress " << NormalStress
697  << ",\n TangentialStress " << TangentialStress
698  << ",\n NormalTraction " << NormalTraction
699  << ",\n TangentialTraction " << TangentialTraction
700  << ",\n Fabric " << Fabric
701  << ",\n CollisionalHeatFlux " << CollisionalHeatFlux
702  << ",\n Dissipation " << Dissipation
703  << ",\n Potential " << Potential
704  << ",\n LocalAngularMomentum " << LocalAngularMomentum
705  << ",\n LocalAngularMomentumFlux " << LocalAngularMomentumFlux
706  << ",\n ContactCoupleStress " << ContactCoupleStress;
707  return ss.str();
708 }
709 
710 template <StatType T>
711 std::string StatisticsPoint<T>::print_sqrt() const
712 {
713  std::stringstream ss;
714  ss<< "Nu " << sqrt(Nu)
715  << ", Density " << sqrt(Density)
716  << ", Momentum " << sqrt(Momentum)
717  << ", DisplacementMomentum " << sqrt(DisplacementMomentum)
718  << ", Displacement " << sqrt(Displacement)
719  << ", MomentumFlux " << sqrt(MomentumFlux)
720  << ", DisplacementMomentumFlux " << sqrt(DisplacementMomentumFlux)
721  << ", EnergyFlux " << sqrt(EnergyFlux)
722  << ", NormalStress " << sqrt(NormalStress)
723  << ", TangentialStress " << sqrt(TangentialStress)
724  << ", NormalTraction " << sqrt(NormalTraction)
725  << ", TangentialTraction " << sqrt(TangentialTraction)
726  << ", Fabric " << sqrt(Fabric)
727  << ", CollisionalHeatFlux " << sqrt(CollisionalHeatFlux)
728  << ", Dissipation " << sqrt(Dissipation)
729  << ", Potential " << sqrt(Potential)
730  << ", LocalAngularMomentum " << LocalAngularMomentum
731  << ", LocalAngularMomentumFlux " << LocalAngularMomentumFlux
732  << ", ContactCoupleStress " << ContactCoupleStress;
733 
734  return ss.str();
735 }
736 
737 template <StatType T>
738 std::string StatisticsPoint<T>::write() const
739 {
740  std::stringstream ss;
741  ss.precision(16);
742  ss<< Nu
743  << " " << Density
744  << " " << Momentum
745  << " " << DisplacementMomentum
746  << " " << Displacement
747  << " " << MomentumFlux
748  << " " << DisplacementMomentumFlux
749  << " " << EnergyFlux
750  << " " << NormalStress
751  << " " << TangentialStress
752  << " " << NormalTraction
753  << " " << TangentialTraction
754  << " " << Fabric
755  << " " << CollisionalHeatFlux
756  << " " << Dissipation
757  << " " << Potential
758  << " " << LocalAngularMomentum
759  << " " << LocalAngularMomentumFlux
760  << " " << ContactCoupleStress;
761  return ss.str();
762 }
763 
764 template <StatType T>
765 void StatisticsPoint<T>::set_Gaussian_invvolume(int dim) {
766  CG_invvolume = compute_Gaussian_invvolume(dim);
767 }
768 
769 template <StatType T>
770 double StatisticsPoint<T>::compute_Gaussian_invvolume(int dim) {
771  if (dim==0) { return 1.; }
772 
773  //this is the prefactor 1/V of phi(r)=1/V*exp(-|r|^2/w^2) for dim=1
774  double CG_invvolume = 1. / ( constants::sqrt_2 * constants::sqrt_pi * get_w());
775  //take into account the cutoff radius and the dimension of the problem
776  if (dim==3) {
777  CG_invvolume = cubic(CG_invvolume);
778  //Wolfram alpha: erf(c/(sqrt(2) w))-(sqrt(2/pi) c e^(-c^2/(2 w^2)))/w
779  CG_invvolume /= erf( get_cutoff() / (constants::sqrt_2*get_w()) )
780  - constants::sqrt_2/constants::sqrt_pi * get_cutoff()/get_w() * exp( -get_cutoff2() / (2*get_w2()) );
781  } else if (dim==2) {
782  CG_invvolume = sqr(CG_invvolume);
783  //Wolfram alpha: integrate(x*exp(-x^2/(2w^2)),{x,0,c})/integrate(x*exp(-x^2/(2w^2)),{x,0,inf})=1-e^(-c^2/(2 w^2))
784  CG_invvolume /= 1-exp( -get_cutoff2() / (2*get_w2()) );
785  } else {
786  CG_invvolume /= erf( get_cutoff() / (constants::sqrt_2*get_w()) );
787  }
788  return CG_invvolume;
789 }
790 
791 template <StatType T>
792 void StatisticsPoint<T>::set_Heaviside_invvolume() {
794  CG_invvolume = 1.0/(4.0/3.0*constants::pi*sqrt(get_w2())*get_w2());
795 }
796 
797 template <StatType T>
798 void StatisticsPoint<T>::set_Polynomial_invvolume(int dim) {
799  if (dim==3) CG_invvolume = 1.0/(sqrt(get_w2())*get_w2());
800  else if (dim==2) CG_invvolume = 1.0/get_w2();
801  else CG_invvolume = 1.0/sqrt(get_w2());
802 }
803 
804 template <StatType T> void StatisticsPoint<T>::set_CG_invvolume() {
805  if (get_CG_type()==HeavisideSphere) {
806  set_Heaviside_invvolume();
807  } else if (get_CG_type()==Gaussian) {
808  set_Gaussian_invvolume(nonaveragedDim());
809  } else if (get_CG_type()==Polynomial) {
810  set_Polynomial_invvolume(nonaveragedDim());
811  } else { std::cerr << "error in CG_function" << std::endl; exit(-1); }
812  CG_invvolume /= averagingVolume();
813 }
814 
815 template<> int StatisticsPoint<XYZ>::nonaveragedDim() { return gb->get_dim_particle(); }
816 template<> int StatisticsPoint<YZ>::nonaveragedDim() { return gb->get_dim_particle()-1; }
817 template<> int StatisticsPoint<XZ>::nonaveragedDim() { return gb->get_dim_particle()-1; }
818 template<> int StatisticsPoint<XY>::nonaveragedDim() { return 2; }
819 template<> int StatisticsPoint<X>::nonaveragedDim() { return 1; }
820 template<> int StatisticsPoint<Y>::nonaveragedDim() { return 1; }
821 template<> int StatisticsPoint<Z>::nonaveragedDim() { return gb->get_dim_particle()-2; }
822 template<> int StatisticsPoint<O>::nonaveragedDim() { return 0; }
823 template<> int StatisticsPoint<RAZ>::nonaveragedDim() { exit(-1); }
824 template<> int StatisticsPoint<AZ>::nonaveragedDim() { exit(-1); }
825 template<> int StatisticsPoint<RZ>::nonaveragedDim() { exit(-1); }
826 template<> int StatisticsPoint<RA>::nonaveragedDim() { exit(-1); }
827 template<> int StatisticsPoint<R>::nonaveragedDim() { exit(-1); }
828 template<> int StatisticsPoint<A>::nonaveragedDim() { exit(-1); }
829 
830 
831 template<> double StatisticsPoint<XYZ>::averagingVolume() {
832  return 1.0;
833 }
834 template<> double StatisticsPoint<XY>::averagingVolume() {
835  return ((gb->get_dim_particle()!=3)?(1.0):(get_zmaxStat()-get_zminStat()));
836 }
837 template<> double StatisticsPoint<XZ>::averagingVolume() {
838  return (get_ymaxStat()-get_yminStat());
839 }
840 template<> double StatisticsPoint<YZ>::averagingVolume() {
841  return (get_xmaxStat()-get_xminStat());
842 }
843 template<> double StatisticsPoint<X>::averagingVolume() {
844  return (get_ymaxStat()-get_yminStat())
845  * ((gb->get_dim_particle()!=3)?(1.0):(get_zmaxStat()-get_zminStat()));
846 }
847 template<> double StatisticsPoint<Y>::averagingVolume() {
848  return (get_xmaxStat()-get_xminStat())
849  * ((gb->get_dim_particle()!=3)?(1.0):(get_zmaxStat()-get_zminStat()));
850 }
851 template<> double StatisticsPoint<Z>::averagingVolume() {
852  return (get_xmaxStat()-get_xminStat())
853  * (get_ymaxStat()-get_yminStat());
854 }
855 template<> double StatisticsPoint<O>::averagingVolume() {
856  return (get_xmaxStat()-get_xminStat())
857  * (get_ymaxStat()-get_yminStat())
858  * ((gb->get_dim_particle()!=3)?(1.0):(get_zmaxStat()-get_zminStat()));
859 }
860 template<> double StatisticsPoint<RAZ>::averagingVolume() { exit(-1); }
861 template<> double StatisticsPoint<AZ>::averagingVolume() { exit(-1); }
862 template<> double StatisticsPoint<RZ>::averagingVolume() { exit(-1); }
863 template<> double StatisticsPoint<RA>::averagingVolume() { exit(-1); }
864 template<> double StatisticsPoint<R>::averagingVolume() { exit(-1); }
865 template<> double StatisticsPoint<A>::averagingVolume() { exit(-1); }
866 
867 template<> Mdouble StatisticsPoint<RA>::CG_function(const Vec3D &PI) {
868  return CG_function_2D(PI);
869 }
870 template<> Mdouble StatisticsPoint<RZ>::CG_function(const Vec3D &PI) {
871  Mdouble dist2 = get_distance2(PI);
872  Mdouble R2 = sqr(PI.X)+sqr(PI.Y);
873  Mdouble RI2 = sqr(Position.X)+sqr(Position.Y);
874  Mdouble Z2 = sqr(PI.Z-Position.Z);
875  if (get_CG_type()==Gaussian) {
876  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume()
877  * exp(-(R2+RI2+Z2)/(2.0*get_w2())) * besselFunc::I0(sqrt(R2*RI2)/get_w2());
878  } else { std::cerr << "error in CG_function<RZ>" << std::endl; exit(-1); }
879 }
880 template<> Mdouble StatisticsPoint<AZ>::CG_function(const Vec3D&) {
881  std::cerr << "error in CG_function<AZ>" << std::endl; exit(-1);
882 }
883 template<> Mdouble StatisticsPoint<R>::CG_function(const Vec3D &PI) {
884  Mdouble dist2 = get_distance2(PI);
885  Mdouble R2 = sqr(PI.X)+sqr(PI.Y);
886  Mdouble RI2 = sqr(Position.X)+sqr(Position.Y);
887  if (get_CG_type()==Gaussian) {
888  return (get_cutoff2()<dist2)?0.0:get_CG_invvolume()
889  * exp(-(R2+RI2)/(2.0*get_w2())) * besselFunc::I0(sqrt(R2*RI2)/get_w2());
890  } else { std::cerr << "error in CG_function<R>" << std::endl; exit(-1); }
891 }
892 template<> Mdouble StatisticsPoint<A>::CG_function(const Vec3D&) {
893  std::cerr << "error in CG_function<A>" << std::endl; exit(-1);
894 }
895 
896 template<> Mdouble StatisticsPoint<XY>::CG_function(const Vec3D &PI) {
897  return CG_function_2D(PI);
898 }
899 template<> Mdouble StatisticsPoint<XZ>::CG_function(const Vec3D &PI) {
900  return CG_function_2D(PI);
901 }
902 template<> Mdouble StatisticsPoint<YZ>::CG_function(const Vec3D &PI) {
903  return CG_function_2D(PI);
904 }
905 template<> Mdouble StatisticsPoint<X>::CG_function(const Vec3D &PI) {
906  return CG_function_1D(PI);
907 }
908 template<> Mdouble StatisticsPoint<Y>::CG_function(const Vec3D &PI) {
909  return CG_function_1D(PI);
910 }
911 template<> Mdouble StatisticsPoint<Z>::CG_function(const Vec3D &PI) {
912  return CG_function_1D(PI);
913 }
914 template<> Mdouble StatisticsPoint<O>::CG_function(const Vec3D&) {
915  return get_CG_invvolume();
916 }
917 
918 template<> Vec3D StatisticsPoint<YZ>::CG_gradient(const Vec3D &P, const Mdouble phi) {
919  //CG_type_todo
920  if (get_CG_type()==Gaussian) {
921  return Vec3D(0.0,P.Y - Position.Y,P.Z - Position.Z)*(phi/get_w2());
922  } else { std::cerr << "error in CG_gradient<YZ>" << std::endl; exit(-1); }
923 }
924 template<> Vec3D StatisticsPoint<XZ>::CG_gradient(const Vec3D &P, const Mdouble phi) {
925  //CG_type_todo
926  if (get_CG_type()==Gaussian) {
927  return Vec3D(P.X - Position.X,0.0,P.Z - Position.Z)*(phi/get_w2());
928  } else { std::cerr << "error in CG_gradient<XZ>" << std::endl; exit(-1); }
929 }
930 template<> Vec3D StatisticsPoint<XY>::CG_gradient(const Vec3D &P, const Mdouble phi) {
931  //CG_type_todo
932  if (get_CG_type()==Gaussian) {
933  return Vec3D(P.X - Position.X,P.Y - Position.Y,0.0)*(phi/get_w2());
934  } else { std::cerr << "error in CG_gradient<XY>" << std::endl; exit(-1); }
935 }
936 
937 template<> Vec3D StatisticsPoint<X>::CG_gradient(const Vec3D &P,const Mdouble phi) {
938  //CG_type_todo
939  if (get_CG_type()==Gaussian) {
940  return Vec3D((P.X - Position.X)*(phi/get_w2()),0.0,0.0);
941  } else if (get_CG_type()==Polynomial) {
942  Mdouble r=fabs(Position.X-P.X)/get_w();
943  return Vec3D(get_CG_invvolume()*evaluatePolynomialGradient(r)/r/get_w2()*(Position.X-P.X),0,0);
944  } else { std::cerr << "error in CG_gradient<X>" << std::endl; exit(-1); }
945 }
946 template<> Vec3D StatisticsPoint<Y>::CG_gradient(const Vec3D &P, const Mdouble phi) {
947  //CG_type_todo
948  if (get_CG_type()==Gaussian) {
949  return Vec3D(0.0,(P.Y - Position.Y)*(phi/get_w2()),0.0);
950  } else if (get_CG_type()==Polynomial) {
951  Mdouble r=fabs(Position.Y-P.Y)/get_w();
952  return Vec3D(0,get_CG_invvolume()*evaluatePolynomialGradient(r)/r/get_w2()*(Position.Y-P.Y),0);
953  } else { std::cerr << "error in CG_gradient<Y>" << std::endl; exit(-1); }
954 }
955 template<> Vec3D StatisticsPoint<Z>::CG_gradient(const Vec3D &P, const Mdouble phi) {
956  //CG_type_todo
957  if (get_CG_type()==Gaussian) {
958  return Vec3D(0.0,0.0,(P.Z - Position.Z)*(phi/get_w2()));
959  } else if (get_CG_type()==Polynomial) {
960  Mdouble r=fabs(Position.Z-P.Z)/get_w();
961  return Vec3D(0,0,get_CG_invvolume()*evaluatePolynomialGradient(r)/r/get_w2()*(Position.Z-P.Z));
962  } else { std::cerr << "error in CG_gradient<Z>" << std::endl; exit(-1); }
963 }
964 
965 template<> Vec3D StatisticsPoint<O>::CG_gradient(const Vec3D&, const Mdouble) {
966  return Vec3D(0.0,0.0,0.0);
967 }
968 
969 //~ template<> Vec3D StatisticsPoint<Z>::CG_integral_gradient(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance) {
970  //~ Vec3D P_P1 = Position - P1;
971  //~ Vec3D P_P2 = Position - P2;
972  //~ Mdouble a = dot(P_P1,P1_P2_normal);
973  //~ Mdouble b = dot(P_P2,P1_P2_normal);
974  //~ Vec3D tangential = P_P1 - a * P1_P2_normal;
975  //~ if (get_CG_type()==Gaussian) {
976  //~ Mdouble wsq2 = sqrt(2*get_w2());
977  //~ Mdouble f = sqrt(2*constants::pi*get_w2());
978  //~ return Vec3D(0.0,0.0,(exp(-sqr(a/wsq2))-exp(-sqr(b/wsq2)))/f/(P2.Z-P1.Z));
979  //~ } else { std::cerr << "error in CG_function" << std::endl; exit(-1); }
980  //~ }
981 
982 template<> Mdouble StatisticsPoint<RA>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi) {
983  Mdouble rpsi_scalar;
984  Mdouble psi = CG_integral_2D(P1, P2, P1_P2_normal, P1_P2_distance, rpsi_scalar);
985  rpsi=Position*psi;
986  return psi;
987 }
988 template<> Mdouble StatisticsPoint<RZ>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D&, Mdouble, Vec3D& rpsi UNUSED) {
989  Vec3D P=(P1+P2)/2; return CG_function(P);
990 }
991 template<> Mdouble StatisticsPoint<AZ>::CG_integral(Vec3D&, Vec3D&, Vec3D&, Mdouble, Vec3D& rpsi UNUSED) {
992  std::cerr << "error in CG_function<AZ>" << std::endl; exit(-1);
993 }
994 template<> Mdouble StatisticsPoint<R>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D&, Mdouble, Vec3D& rpsi UNUSED) {
995  Vec3D P=(P1+P2)/2; return CG_function(P);
996 }
997 template<> Mdouble StatisticsPoint<A>::CG_integral(Vec3D&, Vec3D&, Vec3D&, Mdouble, Vec3D& rpsi UNUSED) {
998  std::cerr << "error in CG_function<A>" << std::endl; exit(-1);
999 }
1000 
1001 template<> Mdouble StatisticsPoint<XY>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi) {
1002  Mdouble rpsi_scalar=0;
1003  Mdouble psi = CG_integral_2D(P1, P2, P1_P2_normal, P1_P2_distance, rpsi_scalar);
1004  rpsi=Vec3D(Position.X*psi,Position.Y*psi,P1.Z*psi-(P1.Z-P2.Z)*rpsi_scalar);
1005  return psi;
1006 }
1007 template<> Mdouble StatisticsPoint<XZ>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi) {
1008  Mdouble rpsi_scalar=0;
1009  Mdouble psi = CG_integral_2D(P1, P2, P1_P2_normal, P1_P2_distance, rpsi_scalar);
1010  rpsi=Vec3D(Position.X*psi,P1.Y*psi-(P1.Y-P2.Y)*rpsi_scalar,Position.Z*psi);
1011  return psi;
1012 }
1013 template<> Mdouble StatisticsPoint<YZ>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi) {
1014  Mdouble rpsi_scalar=0;
1015  Mdouble psi = CG_integral_2D(P1, P2, P1_P2_normal, P1_P2_distance, rpsi_scalar);
1016  rpsi=Vec3D(P1.X*psi-(P1.X-P2.X)*rpsi_scalar,Position.Y*psi,Position.Z*psi);
1017  return psi;
1018 }
1019 template<> Mdouble StatisticsPoint<X>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi) {
1020  Mdouble rpsi_scalar=0;
1021  Mdouble psi = CG_integral_1D(P1, P2, P1_P2_normal, P1_P2_distance, rpsi_scalar);
1022  rpsi=Vec3D(Position.X*psi,P1.Y*psi-(P1.Y-P2.Y)*rpsi_scalar,P1.Z*psi-(P1.Z-P2.Z)*rpsi_scalar);
1023  return psi;
1024 }
1025 template<> Mdouble StatisticsPoint<Y>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi) {
1026  Mdouble rpsi_scalar=0;
1027  Mdouble psi = CG_integral_1D(P1, P2, P1_P2_normal, P1_P2_distance, rpsi_scalar);
1028  rpsi=Vec3D(P1.X*psi-(P1.X-P2.X)*rpsi_scalar,Position.Y*psi,P1.Z*psi-(P1.Z-P2.Z)*rpsi_scalar);
1029  return psi;
1030 }
1031 template<> Mdouble StatisticsPoint<Z>::CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D& rpsi) {
1032  Mdouble rpsi_scalar=0;
1033  Mdouble psi = CG_integral_1D(P1, P2, P1_P2_normal, P1_P2_distance, rpsi_scalar);
1034  rpsi=Vec3D(P1.X*psi-(P1.X-P2.X)*rpsi_scalar,P1.Y*psi-(P1.Y-P2.Y)*rpsi_scalar,Position.Z*psi);
1035  return psi;
1036 }
1037 template<> Mdouble StatisticsPoint<O>::CG_integral(Vec3D &P1 UNUSED, Vec3D &P2 UNUSED, Vec3D &P1_P2_normal UNUSED, Mdouble P1_P2_distance UNUSED, Vec3D& rpsi) {
1038  Mdouble psi = get_CG_invvolume();
1039  rpsi=P1*psi;
1040  return psi;
1041 }
1042 
1043 template<> Vec3D StatisticsPoint<Z>::CG_integral_gradient(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance) {
1044  return Vec3D(0,0,CG_integral_gradient_1D(P1, P2, P1_P2_normal, P1_P2_distance));
1045 }
1046 
1048 template <StatType T> std::ostream& operator<< (std::ostream& os, const StatisticsPoint<T> &stat) {
1049  os.precision(16);
1050  if (stat.mirrorParticle<0) {
1051  //only write std particles, not mirrored particles
1052  os << stat.get_Position() << " " << stat.write() << std::endl;
1053  }
1054  return os;
1055 }
1056 
1057 template <> std::ostream& operator<< (std::ostream& os, const StatisticsPoint<RAZ> &stat) {
1058  os.precision(16); if (stat.mirrorParticle<0) os << stat.get_Position().GetCylindricalCoordinates() << " " << stat.write() << std::endl; return os;
1059 }
1060 template <> std::ostream& operator<< (std::ostream& os, const StatisticsPoint<RA> &stat) {
1061  os.precision(16); if (stat.mirrorParticle<0) os << stat.get_Position().GetCylindricalCoordinates() << " " << stat.write() << std::endl; return os;
1062 }
1063 template <> std::ostream& operator<< (std::ostream& os, const StatisticsPoint<RZ> &stat) {
1064  os.precision(16); if (stat.mirrorParticle<0) os << stat.get_Position().GetCylindricalCoordinates() << " " << stat.write() << std::endl; return os;
1065 }
1066 template <> std::ostream& operator<< (std::ostream& os, const StatisticsPoint<AZ> &stat) {
1067  os.precision(16); if (stat.mirrorParticle<0) os << stat.get_Position().GetCylindricalCoordinates() << " " << stat.write() << std::endl; return os;
1068 }
1069 template <> std::ostream& operator<< (std::ostream& os, const StatisticsPoint<R> &stat) {
1070  os.precision(16); if (stat.mirrorParticle<0) os << stat.get_Position().GetCylindricalCoordinates() << " " << stat.write() << std::endl; return os;
1071 }
1072 template <> std::ostream& operator<< (std::ostream& os, const StatisticsPoint<A> &stat) {
1073  os.precision(16); if (stat.mirrorParticle<0) os << stat.get_Position().GetCylindricalCoordinates() << " " << stat.write() << std::endl; return os;
1074 }
1075 
1076 template<> Mdouble StatisticsPoint<XY>::get_distance2(const Vec3D &P) {return sqr(P.X-Position.X)+sqr(P.Y-Position.Y);}
1077 template<> Mdouble StatisticsPoint<XZ>::get_distance2(const Vec3D &P) {return sqr(P.X-Position.X)+sqr(P.Z-Position.Z);}
1078 template<> Mdouble StatisticsPoint<YZ>::get_distance2(const Vec3D &P) {return sqr(P.Y-Position.Y)+sqr(P.Z-Position.Z);}
1079 template<> Mdouble StatisticsPoint<X> ::get_distance2(const Vec3D &P) {return sqr(P.X-Position.X);}
1080 template<> Mdouble StatisticsPoint<Y> ::get_distance2(const Vec3D &P) {return sqr(P.Y-Position.Y);}
1081 template<> Mdouble StatisticsPoint<Z> ::get_distance2(const Vec3D &P) {return sqr(P.Z-Position.Z);}
1082 template<> Mdouble StatisticsPoint<O> ::get_distance2(const Vec3D&) {return 0;}
1083 
1084 template<> Mdouble StatisticsPoint<XY>::dot(const Vec3D &P, const Vec3D &Q) {return P.X*Q.X + P.Y*Q.Y;}
1085 template<> Mdouble StatisticsPoint<XZ>::dot(const Vec3D &P, const Vec3D &Q) {return P.X*Q.X + P.Z*Q.Z;}
1086 template<> Mdouble StatisticsPoint<YZ>::dot(const Vec3D &P, const Vec3D &Q) {return P.Y*Q.Y + P.Z*Q.Z;}
1087 template<> Mdouble StatisticsPoint<X> ::dot(const Vec3D &P, const Vec3D &Q) {return P.X*Q.X;}
1088 template<> Mdouble StatisticsPoint<Y> ::dot(const Vec3D &P, const Vec3D &Q) {return P.Y*Q.Y;}
1089 template<> Mdouble StatisticsPoint<Z> ::dot(const Vec3D &P, const Vec3D &Q) {return P.Z*Q.Z;}
1090 template<> Mdouble StatisticsPoint<O> ::dot(const Vec3D&, const Vec3D&) {return 0;}
1091 
1092 template<StatType T> Vec3D StatisticsPoint<T>::clearAveragedDirections(Vec3D P) {return P;}
1093 template<> Vec3D StatisticsPoint<XY>::clearAveragedDirections(Vec3D P) {return Vec3D(P.X, P.Y, 0);}
1094 template<> Vec3D StatisticsPoint<XZ>::clearAveragedDirections(Vec3D P) {return Vec3D(P.X, 0, P.Z);}
1095 template<> Vec3D StatisticsPoint<YZ>::clearAveragedDirections(Vec3D P) {return Vec3D(0, P.Y, P.Z);}
1096 template<> Vec3D StatisticsPoint<X>::clearAveragedDirections(Vec3D P) {return Vec3D(P.X, 0, 0);}
1097 template<> Vec3D StatisticsPoint<Y>::clearAveragedDirections(Vec3D P) {return Vec3D(0, P.Y, 0);}
1098 template<> Vec3D StatisticsPoint<Z>::clearAveragedDirections(Vec3D P) {return Vec3D(0, 0, P.Z);}
1099 template<> Vec3D StatisticsPoint<O>::clearAveragedDirections(Vec3D P UNUSED) {return Vec3D(0, 0, 0);}
1100 
1101 template<> Vec3D StatisticsPoint<XY>::cross(Vec3D P, Vec3D &Q) {return Vec3D(0, 0, P.X*Q.Y-P.Y*Q.X);}
1102 template<> Vec3D StatisticsPoint<XZ>::cross(Vec3D P, Vec3D &Q) {return Vec3D(0, P.Z*Q.X-P.X*Q.Z, 0);}
1103 template<> Vec3D StatisticsPoint<YZ>::cross(Vec3D P, Vec3D &Q) {return Vec3D(P.Y*Q.Z-P.Z*Q.Y, 0, 0);}
1104 template<> Vec3D StatisticsPoint<X>::cross(Vec3D P UNUSED, Vec3D &Q UNUSED) {return Vec3D(0,0,0);}
1105 template<> Vec3D StatisticsPoint<Y>::cross(Vec3D P UNUSED, Vec3D &Q UNUSED) {return Vec3D(0,0,0);}
1106 template<> Vec3D StatisticsPoint<Z>::cross(Vec3D P UNUSED, Vec3D &Q UNUSED) {return Vec3D(0,0,0);}
1107 template<> Vec3D StatisticsPoint<O>::cross(Vec3D P UNUSED, Vec3D &Q UNUSED) {return Vec3D(0,0,0);}
1108 
1109 template<> Matrix3D StatisticsPoint<XY>::MatrixCross(Vec3D P, Matrix3D &Q) {return Matrix3D(
1110  0, 0, P.X*Q.YX-P.Y*Q.XX,
1111  0, 0, P.X*Q.YY-P.Y*Q.XY,
1112  0, 0, P.X*Q.YZ-P.Y*Q.XZ);}
1113 template<> Matrix3D StatisticsPoint<XZ>::MatrixCross(Vec3D P, Matrix3D &Q) {return Matrix3D(
1114  0, P.Z*Q.XX-P.X*Q.ZX, 0,
1115  0, P.Z*Q.XY-P.X*Q.ZY, 0,
1116  0, P.Z*Q.XZ-P.X*Q.ZZ, 0);}
1117 template<> Matrix3D StatisticsPoint<YZ>::MatrixCross(Vec3D P, Matrix3D &Q) {return Matrix3D(
1118  P.Y*Q.ZX-P.Z*Q.YX, 0, 0,
1119  P.Y*Q.ZY-P.Z*Q.YY, 0, 0,
1120  P.Y*Q.ZZ-P.Z*Q.YZ, 0, 0);}
1121 template<> Matrix3D StatisticsPoint<X>::MatrixCross(Vec3D P UNUSED, Matrix3D &Q UNUSED) {return Matrix3D(0,0,0,0,0,0,0,0,0);}
1122 template<> Matrix3D StatisticsPoint<Y>::MatrixCross(Vec3D P UNUSED, Matrix3D &Q UNUSED) {return Matrix3D(0,0,0,0,0,0,0,0,0);}
1123 template<> Matrix3D StatisticsPoint<Z>::MatrixCross(Vec3D P UNUSED, Matrix3D &Q UNUSED) {return Matrix3D(0,0,0,0,0,0,0,0,0);}
1124 template<> Matrix3D StatisticsPoint<O>::MatrixCross(Vec3D P UNUSED, Matrix3D &Q UNUSED) {return Matrix3D(0,0,0,0,0,0,0,0,0);}
1125 
Matrix3D::XY
Mdouble XY
Definition: Matrix.h:39
StatisticsPoint::CG_function
Mdouble CG_function(const Vec3D &PI)
Returns the value of the course graining function phi(P,PI)
Definition: StatisticsPoint.hcc:238
StatisticsPoint::get_distance2
Mdouble get_distance2(const Vec3D &P)
returns the coarse graining distance in the coordinates that are not averaged about ...
Definition: StatisticsPoint.hcc:216
StatisticsPoint::DisplacementMomentumFlux
MatrixSymmetric3D DisplacementMomentumFlux
Momentum flux from linear displacement, .
Definition: StatisticsPoint.h:210
Vec3D::X
Mdouble X
Definition: Vector.h:44
Matrix3D::ZY
Mdouble ZY
Definition: Matrix.h:39
StatisticsPoint::ContactCoupleStress
Matrix3D ContactCoupleStress
Definition: StatisticsPoint.h:235
StatisticsPoint::Dissipation
Mdouble Dissipation
Dissipation form collisions, .
Definition: StatisticsPoint.h:228
StatisticsPoint::CG_function_1D
Mdouble CG_function_1D(const Vec3D &PI)
Returns the value of the course graining function phi(P,PI) averaged along a plane.
Definition: StatisticsPoint.hcc:269
Matrix3D::XZ
Mdouble XZ
Definition: Matrix.h:39
Matrix3D::ZX
Mdouble ZX
Definition: Matrix.h:39
constants::sqrt_pi
const Mdouble sqrt_pi
Definition: ExtendedMath.h:55
sqr
#define sqr(a)
Definition: ExtendedMath.h:36
StatisticsPoint::CG_integral_gradient
Vec3D CG_integral_gradient(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance)
gradient of phi
Definition: StatisticsPoint.hcc:293
StatisticsPoint::CG_function_2D
Mdouble CG_function_2D(const Vec3D &PI)
returns the value of the course graining function phi(P,PI) averaged along a line ...
Definition: StatisticsPoint.hcc:251
StatisticsPoint::firstTimeAverage
void firstTimeAverage(const int n)
Defines a division operator needed to time-average values (because the displacement does not have a v...
Definition: StatisticsPoint.hcc:187
operator<<
std::ostream & operator<<(std::ostream &os, const CG S)
Definition: StatisticsPoint.hcc:26
StatisticsPoint::set_Gaussian_invvolume
void set_Gaussian_invvolume(int dim)
sets CG_invvolume if CG_type=Gaussian
Definition: StatisticsPoint.hcc:765
StatisticsPoint::set_CG_invvolume
void set_CG_invvolume()
sets CG_invvolume
Definition: StatisticsPoint.hcc:804
StatisticsPoint::NormalStress
Matrix3D NormalStress
Stress from normal forces, .
Definition: StatisticsPoint.h:214
cubic
#define cubic(a)
Definition: ExtendedMath.h:37
Matrix3D::ZZ
Mdouble ZZ
Definition: Matrix.h:39
Matrix3D::YZ
Mdouble YZ
Definition: Matrix.h:39
StatisticsPoint::set_Heaviside_invvolume
void set_Heaviside_invvolume()
sets CG_invvolume if CG_type=HeaviSideSphere
Definition: StatisticsPoint.hcc:792
StatisticsPoint::CollisionalHeatFlux
Vec3D CollisionalHeatFlux
Heat flux from collisions, .
Definition: StatisticsPoint.h:225
StatisticsPoint::write_variable_names
std::string write_variable_names()
Outputs names of statistical variables in computer-readable format.
Definition: StatisticsPoint.hcc:594
CG
CG
enum used to store the type of coarse-graining function used
Definition: StatisticsPoint.h:27
StatisticsPoint::Nu
Mdouble Nu
Particle volume fraction, .
Definition: StatisticsPoint.h:198
StatisticsPoint::NormalTraction
Vec3D NormalTraction
Traction from normal forces, .
Definition: StatisticsPoint.h:218
StatisticsPoint::CG_integral_2D
Mdouble CG_integral_2D(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Mdouble &rpsi_scalar)
Returns the value of the coarse graining integral averaged along a line.
Definition: StatisticsPoint.hcc:453
StatisticsPoint::DisplacementMomentum
Vec3D DisplacementMomentum
Momentum from linear displacement, , where , with the time intervall between outputs.
Definition: StatisticsPoint.h:204
StatisticsPoint::Fabric
MatrixSymmetric3D Fabric
Fabric tensor, .
Definition: StatisticsPoint.h:222
StatisticsPoint::print
std::string print() const
Outputs statistical variables in human-readable format.
Definition: StatisticsPoint.hcc:685
StatisticsPoint::cross
Vec3D cross(Vec3D P, Vec3D &Q)
Returns the cross product of two vectors in the coordinates that are not averaged about...
Definition: StatisticsPoint.hcc:221
StatisticsPoint::CG_gradient
Vec3D CG_gradient(const Vec3D &P, const Mdouble phi)
gradient of phi
Definition: StatisticsPoint.hcc:281
StatisticsPoint::averagingVolume
double averagingVolume()
StatisticsPoint::EnergyFlux
Vec3D EnergyFlux
Energy flux, .
Definition: StatisticsPoint.h:212
StatisticsPoint::Density
Mdouble Density
Density, .
Definition: StatisticsPoint.h:200
constants::pi
const Mdouble pi
Definition: ExtendedMath.h:54
StatisticsPoint::Displacement
MatrixSymmetric3D Displacement
Linear displacement, .
Definition: StatisticsPoint.h:206
StatisticsPoint::write
std::string write() const
Outputs statistical variables in computer-readable format.
Definition: StatisticsPoint.hcc:738
Mdouble
double Mdouble
Definition: ExtendedMath.h:33
StatisticsPoint::set_zero
void set_zero()
Sets all statistical variables to zero.
Definition: StatisticsPoint.hcc:37
Matrix3D::YX
Mdouble YX
Definition: Matrix.h:39
StatisticsPoint::operator+=
StatisticsPoint< T > & operator+=(const StatisticsPoint< T > &P)
Defines a plus operator needed to average values ( )
Definition: StatisticsPoint.hcc:111
constants::sqrt_2
const Mdouble sqrt_2
Definition: ExtendedMath.h:57
StatisticsPoint::operator/=
StatisticsPoint< T > & operator/=(const Mdouble a)
Defines a division operator needed to average values ( )
Definition: StatisticsPoint.hcc:161
StatisticsPoint::TangentialTraction
Vec3D TangentialTraction
Traction from tangential forces, .
Definition: StatisticsPoint.h:220
StatisticsPoint::LocalAngularMomentumFlux
Matrix3D LocalAngularMomentumFlux
Definition: StatisticsPoint.h:234
Vec3D::Y
Mdouble Y
Definition: Vector.h:44
StatisticsPoint::operator=
StatisticsPoint< T > & operator=(const StatisticsPoint< T > &P)
Defines a equal operator needed for copy constructor.
Definition: StatisticsPoint.hcc:86
StatisticsPoint::clearAveragedDirections
Vec3D clearAveragedDirections(Vec3D P)
Returns a vector where the averaged directions are zero.
Definition: StatisticsPoint.hcc:1092
UNUSED
#define UNUSED
Definition: ExtendedMath.h:38
StatisticsPoint::Potential
Mdouble Potential
Elastic energy .
Definition: StatisticsPoint.h:231
StatisticsPoint::CG_integral_gradient_1D
Mdouble CG_integral_gradient_1D(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance)
Definition: StatisticsPoint.hcc:343
StatisticsPoint
This class stores statistical values for a given spatial position; to be used in combination with Sta...
Definition: StatisticsPoint.h:33
StatisticsPoint::MomentumFlux
MatrixSymmetric3D MomentumFlux
Momentum flux, .
Definition: StatisticsPoint.h:208
Matrix3D::XX
Mdouble XX
Definition: Matrix.h:39
StatisticsPoint::operator-=
StatisticsPoint< T > & operator-=(const StatisticsPoint< T > &P)
Defines a plus operator needed to calculate variance.
Definition: StatisticsPoint.hcc:136
StatisticsVector
This class is used to extract statistical data from MD simulations.
Definition: StatisticsPoint.h:30
StatisticsPoint::getSquared
StatisticsPoint< T > getSquared()
Squares all statistical variables; needed for variance.
Definition: StatisticsPoint.hcc:61
StatisticsPoint::compute_Gaussian_invvolume
double compute_Gaussian_invvolume(int dim)
computes CG_invvolume if CG_type=Gaussian
Definition: StatisticsPoint.hcc:770
StatisticsPoint::CG_integral
Mdouble CG_integral(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Vec3D &rpsi)
Returns the value of the coarse graining integral .
Definition: StatisticsPoint.hcc:400
Vec3D::GetLength
Mdouble GetLength() const
Definition: Vector.h:248
besselFunc::I0
Mdouble I0(Mdouble x)
Definition: ExtendedMath.cc:197
StatisticsPoint::MatrixCross
Matrix3D MatrixCross(Vec3D P, Matrix3D &Q)
Returns the cross product of two vectors in the coordinates that are not averaged about...
Definition: StatisticsPoint.hcc:225
Matrix3D::YY
Mdouble YY
Definition: Matrix.h:39
Matrix3D
Implementation of a 3D matrix.
Definition: Matrix.h:35
Vec3D
Implementation of a 3D vector (by Vitaliy).
Definition: Vector.h:40
StatisticsPoint::Momentum
Vec3D Momentum
Momentum, .
Definition: StatisticsPoint.h:202
Vec3D::Z
Mdouble Z
Definition: Vector.h:44
StatisticsPoint::CG_integral_1D
Mdouble CG_integral_1D(Vec3D &P1, Vec3D &P2, Vec3D &P1_P2_normal, Mdouble P1_P2_distance, Mdouble &rpsi_scalar)
Returns the value of the coarse graining integral averaged along a plane.
Definition: StatisticsPoint.hcc:524
StatisticsPoint::set_Polynomial_invvolume
void set_Polynomial_invvolume(int dim)
sets CG_invvolume if CG_type=Polynomial
Definition: StatisticsPoint.hcc:798
StatisticsPoint::print_sqrt
std::string print_sqrt() const
Outputs root of statistical variables in human-readable format.
Definition: StatisticsPoint.hcc:711
StatisticsPoint::dot
Mdouble dot(const Vec3D &P, const Vec3D &Q)
Returns the dot product of two vectors in the coordinates that are not averaged about.
Definition: StatisticsPoint.hcc:233
StatisticsPoint::TangentialStress
Matrix3D TangentialStress
Stress from tangential forces, .
Definition: StatisticsPoint.h:216