------------------------------------------------------------ revno: 4150 committer: Anton Gladky <[email protected]> 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/[email protected]/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/[email protected]/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/[email protected]/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/[email protected]/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
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp
