------------------------------------------------------------
revno: 3986
committer: raphm1 <raph_mau...@hotmail.com>
timestamp: Wed 2017-01-04 23:18:49 +0100
message:
ViscoelasticPM : fix resolution problem in find_cn_from_en function
Change the value of the initial perturbation for the resolution of the
equation to obtain the damping from the restitution coefficient. Allows for a
much wider range of resolution in terms of mass, stiffness and restitution
coefficient.
Put a limiter to avoid division by zero and force the loop to print the
warning message if "error" goes to nan (no warning in these cases previously)
modified:
pkg/dem/ViscoelasticPM.cpp
--
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 2015-12-16 06:39:37 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2017-01-04 22:18:49 +0000
@@ -19,7 +19,7 @@
/* ViscElPhys */
ViscElPhys::~ViscElPhys(){}
-Real Ip2_ViscElMat_ViscElMat_ViscElPhys::epsilon = 1.0e-15;
+Real Ip2_ViscElMat_ViscElMat_ViscElPhys::epsilon = 1.0e-8;
/* Ip2_ViscElMat_ViscElMat_ViscElPhys */
void Ip2_ViscElMat_ViscElMat_ViscElPhys::go(const shared_ptr<Material>& b1, const shared_ptr<Material>& b2, const shared_ptr<Interaction>& interaction) {
@@ -291,13 +291,16 @@
Real en_temp=get_en_from_cn(cn,m ,kn);
int i =0;
Real error = 1.0/eps;
- while (error > 1.0e-2){
+ while (error > 1.0e-2 or error!=error){
if(i>15){
+ cn = 0.;
+ en_temp = 1.;
cerr<<"Warning in ViscoelasticPM.cpp : Newton-Raphson algorithm did not converged within 15 iterations for contact between "<<interaction->id1<<" and "<<interaction->id2<<". Continue with values : cn="<<cn<<" en="<<en_temp<<endl;
break;
}
i++;
Real deriv=(get_en_from_cn(cn-eps,m ,kn)-get_en_from_cn(cn+eps,m ,kn))/(-2.*eps);
+ deriv = fabs(deriv)>1e-15?deriv:1e-15;
cn = cn - (en_temp-en)/deriv;
en_temp=get_en_from_cn(cn,m ,kn);
error = fabs(en_temp-en)/en;
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : yade-dev@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
