yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3757
committer: Raphael Maurin <raph_maurin@xxxxxxxxxxx>
timestamp: Wed 2013-12-04 10:56:32 +0100
message:
  Change in the calculation of the normal and tangential stiffness and damping. The change affects only the behavior
  when one of the two parameters is zero, e.g. for two particles with different stiffness k1 and k2, the contact
  stiffness will always be k = k1*k2/(k1+k2). Before it gave the same except when k2 = 0 (respectively k1 = 0), where it
  gave k = k1 (resp. k=k2). This is done to ensure continuity in the behavior when one of the two parameter tend
  to zero.
  Add a function contactParameterCalculation in Ip2_ViscElMat_ViscElMat_ViscElPhys to avoid code duplication.
modified:
  pkg/dem/ViscoelasticPM.cpp
  pkg/dem/ViscoelasticPM.hpp


--
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	2013-08-28 13:31:12 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2013-12-04 09:56:32 +0000
@@ -14,8 +14,14 @@
 /* ViscElPhys */
 ViscElPhys::~ViscElPhys(){}
 
+/* Contact parameter calculation function */
+Real Ip2_ViscElMat_ViscElMat_ViscElPhys::contactParameterCalculation(Real l1,Real l2){
+	return (l1>0 or l2>0)?l1*l2/(l1+l2):0;
+}
+
 /* 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) {
+
 	// no updates of an existing contact 
 	if(interaction->phys) return;
 	ViscElMat* mat1 = static_cast<ViscElMat*>(b1.get());
@@ -36,27 +42,18 @@
 	const Real ks2 = mat2->ks*mass2; const Real cs2 = mat2->cs*mass2;
 		
 	ViscElPhys* phys = new ViscElPhys();
-	
-	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;
-	} 
-	
-  if ((mR1>0) or (mR2>0)) {
+
+	phys->kn = contactParameterCalculation(kn1,kn2);
+	phys->ks = contactParameterCalculation(ks1,ks2);
+	phys->cn = contactParameterCalculation(cn1,cn2);
+	phys->cs = contactParameterCalculation(cs1,cs2);
+
+ 	if ((mR1>0) or (mR2>0)) {
 		phys->mR = 2.0/( ((mR1>0)?1/mR1:0) + ((mR2>0)?1/mR2:0) );
 	} else {
 		phys->mR = 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);
 	
@@ -209,3 +206,4 @@
 		addTorque(id2, c2x.cross(f)-momentResistance,scene);
   }
 }
+

=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp	2013-08-28 12:00:42 +0000
+++ pkg/dem/ViscoelasticPM.hpp	2013-12-04 09:56:32 +0000
@@ -65,6 +65,8 @@
 		virtual void go(const shared_ptr<Material>& b1,
 					const shared_ptr<Material>& b2,
 					const shared_ptr<Interaction>& interaction);
+	private :
+		Real contactParameterCalculation(Real l1,Real l2);
 	YADE_CLASS_BASE_DOC(Ip2_ViscElMat_ViscElMat_ViscElPhys,IPhysFunctor,"Convert 2 instances of :yref:`ViscElMat` to :yref:`ViscElPhys` using the rule of consecutive connection.");
 	FUNCTOR2D(ViscElMat,ViscElMat);