yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 4150
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Tue 2014-08-26 16:53:06 +0200
message:
  Fix ViscElPM one more time.
  
  Thanks again to Dominik Boemer.
  http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg08778.html
modified:
  pkg/dem/ViscoelasticPM.cpp
  scripts/checks-and-tests/checks/checkViscElEng.py


--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription

=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2014-08-25 16:21:43 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2014-08-26 14:53:06 +0000
@@ -188,16 +188,19 @@
 		const Real Tc = (tc) ? (*tc)(mat1->id,mat2->id) : (mat1->tc+mat2->tc)/2.0;
 		const Real En = (en) ? (*en)(mat1->id,mat2->id) : (mat1->en+mat2->en)/2.0;
 		const Real Et = (et) ? (*et)(mat1->id,mat2->id) : (mat1->et+mat2->et)/2.0;
-		
-		kn1 = kn2 = 1/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR;
-		cn1 = cn2 = -2.0 /Tc * log(En)*massR;
-		ks1 = ks2 = 2.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR;
-		
-		// It seems to be an error in [Pournin2001] (22) Eq.4, missing factor 2
-		// Thanks to Dominik Boemer for pointing this out
-		// http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg08741.html
-		cs1 = cs2 = -4.0/7.0 /Tc * log(Et)*massR;
-		
+    
+    // Factor 2 at the end of each expression is necessary, because we calculate
+    // individual kn1, kn2, ks1, ks2 etc., because kn1 = 2*kn, ks1 = 2*ks
+    // http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg08778.html
+    kn1 = kn2 = 1/Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(En),2) )*massR*2;
+    cn1 = cn2 = -2.0 /Tc * log(En)*massR*2;
+    ks1 = ks2 = 2.0/7.0 /Tc/Tc * ( Mathr::PI*Mathr::PI + pow(log(Et),2) )*massR*2;
+    cs1 = cs2 = -4.0/7.0 /Tc * log(Et)*massR*2;
+    //           ^^^
+    // It seems to be an error in [Pournin2001] (22) Eq.4, missing factor 2
+    // Thanks to Dominik Boemer for pointing this out
+    // http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg08741.html
+
 		if (std::abs(cn1) <= Mathr::ZERO_TOLERANCE ) cn1=0;
 		if (std::abs(cn2) <= Mathr::ZERO_TOLERANCE ) cn2=0;
 		if (std::abs(cs1) <= Mathr::ZERO_TOLERANCE ) cs1=0;

=== modified file 'scripts/checks-and-tests/checks/checkViscElEng.py'
--- scripts/checks-and-tests/checks/checkViscElEng.py	2014-08-25 15:42:17 +0000
+++ scripts/checks-and-tests/checks/checkViscElEng.py	2014-08-26 14:53:06 +0000
@@ -45,13 +45,13 @@
 
   if v0<=0 and v>0:
     en=-v/v0
-    print ("Precalculated en value %f" % 0.734839832393159)
-    print ("Obtained en value %.15f" % en)
+    print ("Precalculated en value %.12f" % 0.736356797441)
+    print ("Obtained en value %.12f" % en)
     O.pause()
   v0=v
 
 O.run(1000000)
 O.wait()
 
-if ((abs(0.734839832393159-en)/en)>tolerance):
+if ((abs(0.736356797441-en)/en)>tolerance):
   resultStatus += 1