yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3986
committer: raphm1 <raph_maurin@xxxxxxxxxxx>
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;