yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3872
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Wed 2014-04-02 17:33:41 +0200
message:
  Deprecate getViscoelasticFromSpheresInteraction.
  
  From now all parameter calculations in Viscolastic contact
  model will be executed directly in Ip2_ViscElMat_ViscElMat_ViscElPhys.
  
  Using of external function getViscoelasticFromSpheresInteraction is
  deprecated now and should be avoided. Minor update of scripts,
  which are using it is necessary.
modified:
  pkg/dem/Shop_02.cpp
  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/Shop_02.cpp'
--- pkg/dem/Shop_02.cpp	2014-02-27 08:27:24 +0000
+++ pkg/dem/Shop_02.cpp	2014-04-02 15:33:41 +0000
@@ -147,13 +147,7 @@
 
 void Shop::getViscoelasticFromSpheresInteraction( Real tc, Real en, Real es, shared_ptr<ViscElMat> b)
 {
-    b->kn = 1/tc/tc * ( Mathr::PI*Mathr::PI + pow(log(en),2) );
-    b->cn = -2.0 /tc * log(en);
-    b->ks = 2.0/7.0 /tc/tc * ( Mathr::PI*Mathr::PI + pow(log(es),2) );
-    b->cs = -2.0/7.0 /tc * log(es);
-
-    if (abs(b->cn) <= Mathr::ZERO_TOLERANCE ) b->cn=0;
-    if (abs(b->cs) <= Mathr::ZERO_TOLERANCE ) b->cs=0;
+    throw runtime_error("Setting parameters in ViscoElastic model is changed. You do not need to use getViscoelasticFromSpheresInteraction function any more, because this functino is deprecated. You need to set the parameters tc, en and es directly in material properties. Please, update your scripts.");
 }
 
 /* This function is copied almost verbatim from scientific python, module Visualization, class ColorScale

=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2014-04-01 14:45:46 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2014-04-02 15:33:41 +0000
@@ -124,38 +124,90 @@
 	Real mass1 = 1.0;
 	Real mass2 = 1.0;
 	
-	if (mat1->massMultiply and mat2->massMultiply) {
-		mass1 = Body::byId(interaction->getId1())->state->mass;
-		mass2 = Body::byId(interaction->getId2())->state->mass;
-		if (mass1 == 0.0 and mass2 > 0.0) {
-			mass1 = mass2;
-		} else if (mass2 == 0.0 and mass1 > 0.0) {
-			mass2 = mass1;
-		}
+	if ((isnormal(mat1->kn) and  not (isnormal(mat2->kn))) or
+			(isnormal(mat2->kn) and  not (isnormal(mat1->kn))) or
+			(isnormal(mat1->ks) and  not (isnormal(mat2->ks))) or
+			(isnormal(mat2->ks) and  not (isnormal(mat1->ks))) or
+			(isnormal(mat1->cn) and  not (isnormal(mat2->cn))) or
+			(isnormal(mat2->cn) and  not (isnormal(mat1->cn))) or
+			(isnormal(mat1->cs) and  not (isnormal(mat2->cs))) or
+			(isnormal(mat2->cs) and  not (isnormal(mat1->cs))) or
+			(isnormal(mat1->tc) and  not (isnormal(mat2->tc))) or
+			(isnormal(mat2->tc) and  not (isnormal(mat1->tc))) or
+			(isnormal(mat1->en) and  not (isnormal(mat2->en))) or
+			(isnormal(mat2->en) and  not (isnormal(mat1->en))) or
+			(isnormal(mat1->et) and  not (isnormal(mat2->et))) or
+			(isnormal(mat2->et) and  not (isnormal(mat1->et)))) {
+				throw runtime_error("Both materials should have the same defined set of variables e.g. tc, ks etc.!"); 
+			}
+			
+	mass1 = Body::byId(interaction->getId1())->state->mass;
+	mass2 = Body::byId(interaction->getId2())->state->mass;
+	if (mass1 == 0.0 and mass2 > 0.0) {
+		mass1 = mass2;
+	} else if (mass2 == 0.0 and mass1 > 0.0) {
+		mass2 = mass1;
 	}
 	
+	const Real massR = 2*mass1*mass2/(mass1+mass2);
+	
 	GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
 	Real R1=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
 	Real R2=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
 	
-	const Real kn1 = isnan(mat1->kn)?2*mat1->young*R1:mat1->kn*mass1;
-	const Real cn1 = mat1->cn*mass1;
-	const Real ks1 = isnan(mat1->ks)?kn1*mat1->poisson:mat1->ks*mass1;
-	const Real cs1 = mat1->cs*mass1;
+	Real kn1 = 0.0; Real kn2 = 0.0;
+	Real cn1 = 0.0; Real cn2 = 0.0;
+	Real ks1 = 0.0; Real ks2 = 0.0;
+	Real cs1 = 0.0; Real cs2 = 0.0;
+	
+	if ((isnormal(mat1->tc)) and (isnormal(mat1->en)) and (isnormal(mat1->et))) {
+		//Set parameters according to [Pournin2001]
+		
+		const Real tc = (mat1->tc+mat2->tc)/2.0;
+		const Real en = (mat1->en+mat2->en)/2.0;
+		const Real et = (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;
+		cs1 = cs2 = -2.0/7.0 /tc * log(et)*massR;
+	
+		if (abs(cn1) <= Mathr::ZERO_TOLERANCE ) cn1=0;
+		if (abs(cn2) <= Mathr::ZERO_TOLERANCE ) cn2=0;
+		if (abs(cs1) <= Mathr::ZERO_TOLERANCE ) cs1=0;
+		if (abs(cs2) <= Mathr::ZERO_TOLERANCE ) cs2=0;
+	} else if ((isnormal(mat1->kn)) and (isnormal(mat1->ks)) and (isnormal(mat1->cn)) and (isnormal(mat1->cs))) {
+		//Set parameters explicitly
+		kn1 = mat1->kn;
+		kn2 = mat2->kn;
+		ks1 = mat1->ks;
+		ks2 = mat2->ks;
+		cn1 = mat1->cn;
+		cn2 = mat2->cn;
+		cs1 = mat1->cs;
+		cs2 = mat2->cs;
+	} else {
+		//Set parameters on the base of young modulus
+		kn1 = 2*mat1->young*R1;
+		kn2 = 2*mat2->young*R2;
+		ks1 = kn1*mat1->poisson;
+		ks2 = kn2*mat2->poisson;
+		if ((isnormal(mat1->cn)) and (isnormal(mat1->cs))) {
+			cn1 = mat1->cn;
+			cn2 = mat2->cn;
+			cs1 = mat1->cs;
+			cs2 = mat2->cs;
+		}
+	}
 	
 	const Real mR1 = mat1->mR;      const Real mR2 = mat2->mR; 
 	const int mRtype1 = mat1->mRtype; const int mRtype2 = mat2->mRtype;
 	
-	const Real kn2 = isnan(mat2->kn)?2*mat2->young*R2:mat2->kn*mass2;
-	const Real cn2 = mat2->cn*mass2;
-	const Real ks2 = isnan(mat2->ks)?kn2*mat2->poisson:mat2->ks*mass2;
-	const Real cs2 = mat2->cs*mass2;
 	
-
-	phys->kn = contactParameterCalculation(kn1,kn2, mat1->massMultiply&&mat2->massMultiply);
-	phys->ks = contactParameterCalculation(ks1,ks2, mat1->massMultiply&&mat2->massMultiply);
-	phys->cn = contactParameterCalculation(cn1,cn2, mat1->massMultiply&&mat2->massMultiply);
-	phys->cs = contactParameterCalculation(cs1,cs2, mat1->massMultiply&&mat2->massMultiply);
+	phys->kn = contactParameterCalculation(kn1,kn2, 1);
+	phys->ks = contactParameterCalculation(ks1,ks2, 1);
+	phys->cn = contactParameterCalculation(cn1,cn2, 1);
+	phys->cs = contactParameterCalculation(cs1,cs2, 1);
 
  	if ((mR1>0) or (mR2>0)) {
 		phys->mR = 2.0/( ((mR1>0)?1/mR1:0) + ((mR2>0)?1/mR2:0) );

=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp	2014-04-01 14:45:46 +0000
+++ pkg/dem/ViscoelasticPM.hpp	2014-04-02 15:33:41 +0000
@@ -18,11 +18,13 @@
 	public:
 		virtual ~ViscElMat();
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscElMat,FrictMat,"Material for simple viscoelastic model of contact.\n\n.. note::\n\t ``Shop::getViscoelasticFromSpheresInteraction`` (and :yref:`yade.utils.getViscoelasticFromSpheresInteraction` in python) compute :yref:`kn<ViscElMat::kn>`, :yref:`cn<ViscElMat::cn>`,  :yref:`ks<ViscElMat::ks>`,  :yref:`cs<ViscElMat::cs>` from analytical solution of a pair spheres interaction problem.",
-		((Real,kn,NaN,,"Normal elastic stiffness"))
-		((Real,cn,NaN,,"Normal viscous constant"))
-		((Real,ks,NaN,,"Shear elastic stiffness"))
-		((Real,cs,NaN,,"Shear viscous constant"))
-		((bool,massMultiply,true,,"Stiffness and viscosity are multiplied by the reduced mass. If massMultiply=false, these parameter are set explicitly without mass multiplication"))
+		((Real,tc,NaN,,"Contact time"))
+		((Real,en,NaN,,"Restitution coefficient in normal direction"))
+		((Real,et,NaN,,"Restitution coefficient in tangential direction"))
+		((Real,kn,NaN,,"Normal elastic stiffness. Attention, this parameter cannot be set if tc, en or es is defined!"))
+		((Real,cn,NaN,,"Normal viscous constant. Attention, this parameter cannot be set if tc, en or es is defined!"))
+		((Real,ks,NaN,,"Shear elastic stiffness. Attention, this parameter cannot be set if tc, en or es is defined!"))
+		((Real,cs,NaN,,"Shear viscous constant. Attention, this parameter cannot be set if tc, en or es is defined!"))
 		((Real,mR,0.0,,"Rolling resistance, see [Zhou1999536]_."))
 		((unsigned int,mRtype,1,,"Rolling resistance type, see [Zhou1999536]_. mRtype=1 - equation (3) in [Zhou1999536]_; mRtype=2 - equation (4) in [Zhou1999536]_."))
 		((bool,Capillar,false,,"True, if capillar forces need to be added."))