← Back to team overview

yade-dev team mailing list archive

[Branch ~yade-dev/yade/trunk] Rev 2879: Fixes nans in ViscoelasticPM, when kn=0 or ks=0. Fixes 804650

 

------------------------------------------------------------
revno: 2879
fixes bug(s): https://launchpad.net/bugs/804650
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
branch nick: yade
timestamp: Sat 2011-07-02 18:46:27 +0200
message:
  Fixes nans in ViscoelasticPM, when kn=0 or ks=0. Fixes 804650
modified:
  pkg/dem/ViscoelasticPM.cpp


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

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2011-01-11 14:09:30 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2011-07-02 16:46:27 +0000
@@ -26,10 +26,21 @@
 	const Real kn2 = mat2->kn*mass2; const Real cn2 = mat2->cn*mass2;
 	const Real ks2 = mat2->ks*mass2; const Real cs2 = mat2->cs*mass2;
 	ViscElPhys* phys = new ViscElPhys();
-	phys->kn = 1/( (kn1?1/kn1:0) + (kn2?1/kn2:0) );
-	phys->ks = 1/( (ks1?1/ks1:0) + (ks2?1/ks2:0) );
+	
+	if ((kn1>0) or (kn2>0)) {
+		phys->kn = 1/( ((kn1>0)?1/kn1:0) + ((kn2>0)?1/kn2:0) );
+	} else {
+		phys->kn = 0;
+	}
+	if ((ks1>0) or (ks2>0)) {
+		phys->ks = 1/( ((ks1>0)?1/ks1:0) + ((ks2>0)?1/ks2:0) );
+	} else {
+		phys->ks = 0;
+	} 
+	
 	phys->cn = (cn1?1/cn1:0) + (cn2?1/cn2:0); phys->cn = phys->cn?1/phys->cn:0;
 	phys->cs = (cs1?1/cs1:0) + (cs2?1/cs2:0); phys->cs = phys->cs?1/phys->cs:0;
+	
 	phys->tangensOfFrictionAngle = std::tan(std::min(mat1->frictionAngle, mat2->frictionAngle)); 
 	phys->shearForce = Vector3r(0,0,0);
 	//phys->prevNormal = Vector3r(0,0,0);