yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 3797
committer: Francois <francois.kneib@xxxxxxxxx>
timestamp: Wed 2014-01-08 15:36:59 +0100
message:
  Allow the contact between (coh)frictMat and viscElMat as frictPhys. Almost everything was done by inheritance, just had to convert stiffnesses to modulus and modulus to stiffnesses to ensure material compatibility.
  Note that for the moment the timeStepper cannot handle this kind of simulations -> will be fixed soon.
modified:
  pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.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/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.cpp	2013-12-18 14:45:51 +0000
+++ pkg/dem/Ip2_FrictMat_FrictMat_FrictPhys.cpp	2014-01-08 14:36:59 +0000
@@ -13,6 +13,7 @@
 #include<yade/core/Omega.hpp>
 #include<yade/core/Scene.hpp>
 #include<yade/pkg/common/ElastMat.hpp>
+#include<yade/pkg/dem/ViscoelasticPM.hpp>
 #include <cassert>
 
 
@@ -22,8 +23,32 @@
 					, const shared_ptr<Interaction>& interaction)
 {
 	if(interaction->phys) return;
+	
 	const shared_ptr<FrictMat>& mat1 = YADE_PTR_CAST<FrictMat>(b1);
 	const shared_ptr<FrictMat>& mat2 = YADE_PTR_CAST<FrictMat>(b2);
+	
+	Real Ra,Rb;//Vector3r normal;
+	assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode
+	GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
+	Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
+	Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
+	
+	//The two contitions above are used in the case of a ViscElMat/FrictMat contact, if kn and ks are set for the ViscElMat. Young and Poisson are deduced. This is not a perfect solution and the user should set young/poisson only.
+	if( b1->getClassIndex()==ViscElMat::getClassIndexStatic()){
+		const shared_ptr<ViscElMat>& mat1 = YADE_PTR_CAST<ViscElMat>(b1);
+		if( (!isnan(mat1->kn)) && (!isnan(mat1->ks)) ){
+			mat1->young=mat1->kn/(2.*Ra);
+			mat1->poisson=mat1->ks/mat1->kn;
+		}
+	}
+	if( b2->getClassIndex()==ViscElMat::getClassIndexStatic()){
+		const shared_ptr<ViscElMat>& mat2 = YADE_PTR_CAST<ViscElMat>(b2);
+		if( (!isnan(mat2->kn)) && (!isnan(mat2->ks)) ){
+			mat2->young=mat2->kn/(2.*Rb);
+			mat2->poisson=mat2->ks/mat2->kn;
+		}
+	}
+	
 	interaction->phys = shared_ptr<FrictPhys>(new FrictPhys());
 	const shared_ptr<FrictPhys>& contactPhysics = YADE_PTR_CAST<FrictPhys>(interaction->phys);
 	Real Ea 	= mat1->young;
@@ -31,12 +56,6 @@
 	Real Va 	= mat1->poisson;
 	Real Vb 	= mat2->poisson;
 
-	Real Ra,Rb;//Vector3r normal;
-	assert(dynamic_cast<GenericSpheresContact*>(interaction->geom.get()));//only in debug mode
-	GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->geom.get());
-	Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2;
-	Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1;
-
 	//harmonic average of the two stiffnesses when (2*Ri*Ei) is the stiffness of a contact point on sphere "i"
 	Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
 	//same for shear stiffness

=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp	2014-01-07 15:37:27 +0000
+++ pkg/dem/ViscoelasticPM.cpp	2014-01-08 14:36:59 +0000
@@ -42,12 +42,22 @@
 		mass2 = Body::byId(interaction->getId2())->state->mass;
 	}
 	
-	const Real kn1 = mat1->kn*mass1; const Real cn1 = mat1->cn*mass1;
-	const Real ks1 = mat1->ks*mass1; const Real cs1 = mat1->cs*mass1;
+	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;
+	
 	const Real mR1 = mat1->mR;      const Real mR2 = mat2->mR; 
 	const int mRtype1 = mat1->mRtype; const int mRtype2 = mat2->mRtype;
-	const Real kn2 = mat2->kn*mass2; const Real cn2 = mat2->cn*mass2;
-	const Real ks2 = mat2->ks*mass2; const Real cs2 = mat2->cs*mass2;
+	
+	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;
 		
 	ViscElPhys* phys = new ViscElPhys();
 
@@ -111,6 +121,7 @@
 	
 	if (geom.penetrationDepth<0) {
 		if (phys.liqBridgeCreated and -geom.penetrationDepth<phys.sCrit and phys.Capillar) {
+			cout<<id1<<" "<<id2<<endl;
 			phys.normalForce = -calculateCapillarForce(geom, phys)*geom.normal;
 		  if (I->isActive) {
 				addForce (id1,-phys.normalForce,scene);

=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp	2013-12-18 08:04:27 +0000
+++ pkg/dem/ViscoelasticPM.hpp	2014-01-08 14:36:59 +0000
@@ -4,24 +4,25 @@
 
 #pragma once
 
-#include<yade/core/Material.hpp>
+#include<yade/pkg/common/ElastMat.hpp>
 #include<yade/pkg/dem/FrictPhys.hpp>
 #include<yade/pkg/common/Dispatching.hpp>
 #include<yade/pkg/dem/ScGeom.hpp>
+#include<yade/pkg/dem/DemXDofGeom.hpp>
+
 
 /* Simple viscoelastic model */
 
 /// Material
 /// Note: Shop::getViscoelasticFromSpheresInteraction can get kn,cn,ks,cs from a analytical solution of a pair spheres interaction problem.
-class ViscElMat : public Material {
+class ViscElMat : public FrictMat {
 	public:
 		virtual ~ViscElMat();
-	YADE_CLASS_BASE_DOC_ATTRS_CTOR(ViscElMat,Material,"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.",
+	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"))
-		((Real,frictionAngle,NaN,,"Friction angle [rad]"))
 		((bool,massMultiply,true,,"Stiffness and viscosity are multiplied by the reduced mass. If massMultiply=false, these parameter are set explicitly without mass multiplication"))
 		((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]_."))
@@ -32,7 +33,7 @@
 		((std::string,CapillarType,"",,"Different types of capillar interaction: Willett_numeric, Willett_analytic [Willett2000]_ , Weigert [Weigert1999]_ , Rabinovich [Rabinov2005]_ ")),
 		createIndex();
 	);
-	REGISTER_CLASS_INDEX(ViscElMat,Material);
+	REGISTER_CLASS_INDEX(ViscElMat,FrictMat);
 };
 REGISTER_SERIALIZABLE(ViscElMat);
 
@@ -81,7 +82,9 @@
 		virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
 	private:
 		Real calculateCapillarForce(const ScGeom& geom, ViscElPhys& phys);
+	public :
 	FUNCTOR2D(ScGeom,ViscElPhys);
 	YADE_CLASS_BASE_DOC(Law2_ScGeom_ViscElPhys_Basic,LawFunctor,"Linear viscoelastic model operating on :yref:`ScGeom` and :yref:`ViscElPhys`. The model is mostly based on the paper for For details see Pournin [Pournin2001]_ .");
+	DECLARE_LOGGER;
 };
 REGISTER_SERIALIZABLE(Law2_ScGeom_ViscElPhys_Basic);