yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2509
committer: bchareyre <bchareyre@dt-rv020>
branch nick: yade
timestamp: Mon 2010-10-25 16:52:21 +0200
message:
  -Move cohesion and moment parameters from functors to bodies and interactions, and define them in Ip2_2xCohFrictMat_...
modified:
  pkg/dem/CohFrictMat.hpp
  pkg/dem/CohFrictPhys.hpp
  pkg/dem/CohesiveFrictionalContactLaw.cpp
  pkg/dem/CohesiveFrictionalContactLaw.hpp
  pkg/dem/CohesiveTriaxialTest.cpp
  pkg/dem/ElasticContactLaw.cpp
  pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp
  pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.hpp
  scripts/test/chained-cylinder-spring.py


--
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/CohFrictMat.hpp'
--- pkg/dem/CohFrictMat.hpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/CohFrictMat.hpp	2010-10-25 14:52:21 +0000
@@ -21,7 +21,10 @@
 /// Serialization
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CohFrictMat,FrictMat,"",
 		((bool,isBroken,true,,""))
-		((bool,isCohesive,true,,"")),
+		((bool,isCohesive,true,,""))
+		((Real,normalCohesion,10000000,,""))
+		((Real,shearCohesion,10000000,,""))
+		((bool,momentRotationLaw,false,,"Use bending/twisting moment at contact. The contact will have moments only if both bodies have this flag true. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details.")),
 		createIndex();
 					);
 /// Indexable

=== modified file 'pkg/dem/CohFrictPhys.hpp'
--- pkg/dem/CohFrictPhys.hpp	2010-10-20 11:31:49 +0000
+++ pkg/dem/CohFrictPhys.hpp	2010-10-25 14:52:21 +0000
@@ -24,9 +24,11 @@
 		((Real,kr,0,,"rotational stiffness [N.m/rad]"))
 		((Real,normalAdhesion,0,,"tensile strength"))
 		((Real,shearAdhesion,0,,"cohesive part of the shear strength (a frictional term might be added depending on :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment::always_use_moment_law`)"))
+		((bool,momentRotationLaw,false,,"use bending/twisting moment at contacts. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))
+		((Real,creep_viscosity,-1,,"creep viscosity [Pa.s/m]."))
 		// internal attributes
-		((Vector3r,moment_twist,Vector3r(0,0,0),,""))
-		((Vector3r,moment_bending,Vector3r(0,0,0),,""))
+		((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment"))
+		((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment"))
 		,
 		createIndex();
 	);

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp	2010-10-20 11:31:49 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp	2010-10-25 14:52:21 +0000
@@ -13,28 +13,22 @@
 #include<yade/core/Omega.hpp>
 #include<yade/core/Scene.hpp>
 
-
 YADE_PLUGIN((CohesiveFrictionalContactLaw)(Law2_ScGeom_CohFrictPhys_CohesionMoment));
 CREATE_LOGGER(Law2_ScGeom_CohFrictPhys_CohesionMoment);
 
 Vector3r translation_vect_ ( 0.10,0,0 );
 
-
-
 void CohesiveFrictionalContactLaw::action()
 {
 	if(!functor) functor=shared_ptr<Law2_ScGeom_CohFrictPhys_CohesionMoment>(new Law2_ScGeom_CohFrictPhys_CohesionMoment);
 	functor->always_use_moment_law = always_use_moment_law;
 	functor->shear_creep=shear_creep;
 	functor->twist_creep=twist_creep;
-	functor->creep_viscosity = creep_viscosity;
+	functor->creep_viscosity=creep_viscosity;
 	functor->scene=scene;
-	functor->momentRotationLaw=momentRotationLaw;
-	
 	FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
 		if(!I->isReal()) continue;
-		functor->go(I->geom, I->phys, I.get());
-	}
+		functor->go(I->geom, I->phys, I.get());}
 }
 
 void Law2_ScGeom_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
@@ -53,7 +47,6 @@
 	if (un < 0 && (currentContactPhysics->normalForce.squaredNorm() > pow(currentContactPhysics->normalAdhesion,2)
 	               || currentContactPhysics->normalAdhesion==0)) {
 		// BREAK due to tension
-		
 		scene->interactions->requestErase(contact->getId1(),contact->getId2());
 	} else {
 		State* de1 = Body::byId(id1,scene)->state.get();
@@ -85,7 +78,7 @@
 		applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
 
 		/// Moment law        ///
-		if (momentRotationLaw && (!currentContactPhysics->cohesionBroken || always_use_moment_law)) {
+		if (currentContactPhysics->momentRotationLaw && (!currentContactPhysics->cohesionBroken || always_use_moment_law)) {
 			if (twist_creep) {
 				Real viscosity_twist = creep_viscosity * std::pow((2 * std::min(currentContactGeometry->radius1,currentContactGeometry->radius2)),2) / 16.0;
 				Real angle_twist_creeped = currentContactGeometry->getTwist() * (1 - dt/viscosity_twist);
@@ -101,6 +94,5 @@
 			scene->forces.addTorque(id2, moment);
 		}
 		/// Moment law END       ///
-		currentContactPhysics->prevNormal = currentContactGeometry->normal;
 	}
 }

=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp	2010-10-20 11:31:49 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp	2010-10-25 14:52:21 +0000
@@ -16,9 +16,8 @@
 class Law2_ScGeom_CohFrictPhys_CohesionMoment: public LawFunctor{
 	public:
 	virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
-	YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. There is no maximum value of moments in the current implementation, though they could be added in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
+	YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_n$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. There is no maximum value of moments in the current implementation, though they could be added in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.\n\n.. note::\n  Periodicity is not handled yet in this law.",
 		((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
-		((bool,momentRotationLaw,false,,"use bending/twisting moment at contacts. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))
 		((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
 		((bool,shear_creep,false,,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
 		((bool,twist_creep,false,,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
@@ -40,7 +39,6 @@
 		
 	YADE_CLASS_BASE_DOC_ATTRS(CohesiveFrictionalContactLaw,GlobalEngine,"[DEPRECATED] Loop over interactions applying :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment` on all interactions.\n\n.. note::\n  Use :yref:`InteractionLoop` and :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment` instead of this class for performance reasons.",
 		((bool,neverErase,false,,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
-		((bool,momentRotationLaw,false,,"use bending/twisting moment at contacts. See :yref:`CohesiveFrictionalContactLaw::always_use_moment_law` for details."))
 		((bool,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
 		((bool,shear_creep,false,,"activate creep on the shear force, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))
 		((bool,twist_creep,false,,"activate creep on the twisting moment, using :yref:`CohesiveFrictionalContactLaw::creep_viscosity`."))

=== modified file 'pkg/dem/CohesiveTriaxialTest.cpp'
--- pkg/dem/CohesiveTriaxialTest.cpp	2010-10-20 11:31:49 +0000
+++ pkg/dem/CohesiveTriaxialTest.cpp	2010-10-25 14:52:21 +0000
@@ -224,7 +224,9 @@
 	physics->young			= sphereYoungModulus;
 	physics->poisson		= sphereKsDivKn;
 	physics->frictionAngle		= sphereFrictionDeg * Mathr::PI/180.0;
-
+	physics->shearCohesion = shearCohesion;
+	physics->normalCohesion = normalCohesion;
+	
 	if((!dynamic) && (!boxWalls))
 	{
 		physics->young			= boxYoungModulus;
@@ -270,7 +272,9 @@
 	physics->young			= boxYoungModulus;
 	physics->poisson		= boxKsDivKn;
 	physics->frictionAngle		= boxFrictionDeg * Mathr::PI/180.0;
-
+	physics->shearCohesion = 0;
+	physics->normalCohesion = 0;
+	
 	aabb->color		= Vector3r(1,0,0);
 	
 	iBox->extents			= extents;
@@ -279,7 +283,7 @@
 
 	body->bound		= aabb;
 	body->shape	= iBox;
-	body->material	= physics;	
+	body->material	= physics;
 }
 
 
@@ -293,8 +297,6 @@
 	interactionGeometryDispatcher->add(s2);
 
 	shared_ptr<Ip2_2xCohFrictMat_CohFrictPhys> cohesiveFrictionalRelationships = shared_ptr<Ip2_2xCohFrictMat_CohFrictPhys> (new Ip2_2xCohFrictMat_CohFrictPhys);
-	cohesiveFrictionalRelationships->shearCohesion = shearCohesion;
-	cohesiveFrictionalRelationships->normalCohesion = normalCohesion;
 	cohesiveFrictionalRelationships->setCohesionOnNewContacts = setCohesionOnNewContacts;
 	shared_ptr<IPhysDispatcher> interactionPhysicsDispatcher(new IPhysDispatcher);
 	interactionPhysicsDispatcher->add(cohesiveFrictionalRelationships);

=== modified file 'pkg/dem/ElasticContactLaw.cpp'
--- pkg/dem/ElasticContactLaw.cpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/ElasticContactLaw.cpp	2010-10-25 14:52:21 +0000
@@ -94,8 +94,6 @@
 		scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
 		scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
 	}
-	//FIXME : make sure phys->prevNormal is not used anywhere, remove it from physics and remove this line :
-	phys->prevNormal = geom->normal;
 }
 
 // same as elasticContactLaw, but using Dem3DofGeom

=== modified file 'pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp	2010-10-20 11:31:49 +0000
+++ pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.cpp	2010-10-25 14:52:21 +0000
@@ -57,13 +57,14 @@
 			if ((setCohesionOnNewContacts || setCohesionNow) && sdec1->isCohesive && sdec2->isCohesive)
 			{
 				contactPhysics->cohesionBroken = false;
-				contactPhysics->normalAdhesion			= normalCohesion*pow(std::min(Db, Da),2);
-				contactPhysics->shearAdhesion			= shearCohesion*pow(std::min(Db, Da),2);;
+				contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(Db, Da),2);
+				contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(Db, Da),2);
 				geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
 			}
 			contactPhysics->kn = Kn;
 			contactPhysics->ks = Ks;
 			contactPhysics->kr = Kr;
+			contactPhysics->momentRotationLaw=(sdec1->momentRotationLaw && sdec2->momentRotationLaw);
 			//contactPhysics->elasticRollingLimit = elasticRollingLimit;
 
 		}
@@ -72,8 +73,8 @@
 			if (setCohesionNow && sdec1->isCohesive && sdec2->isCohesive)
 			{
 				contactPhysics->cohesionBroken = false;
-				contactPhysics->normalAdhesion			= normalCohesion*pow(std::min(geom->radius2, geom->radius1),2);
-				contactPhysics->shearAdhesion			= shearCohesion*pow(std::min(geom->radius2, geom->radius1),2);
+				contactPhysics->normalAdhesion = std::min(sdec1->normalCohesion,sdec2->normalCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
+				contactPhysics->shearAdhesion = std::min(sdec1->shearCohesion,sdec2->shearCohesion)*pow(std::min(geom->radius2, geom->radius1),2);
 				geom->initRotations(*(Body::byId(interaction->getId1(),scene)->state),*(Body::byId(interaction->getId2(),scene)->state));
 			}
 		}

=== modified file 'pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.hpp'
--- pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.hpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/Ip2_2xCohFrictMat_CohFrictPhys.hpp	2010-10-25 14:52:21 +0000
@@ -22,8 +22,6 @@
 		"Generates cohesive-frictional interactions with moments. Used in the contact law :yref:`Law2_ScGeom_CohFrictPhys_CohesionMoment`.",
 		((bool,setCohesionNow,false,,""))
 		((bool,setCohesionOnNewContacts,false,,""))
-		((Real,normalCohesion,10000000,,""))
-		((Real,shearCohesion,10000000,,""))
 		,
 		cohesionDefinitionIteration = -1; 
 		);

=== modified file 'scripts/test/chained-cylinder-spring.py'
--- scripts/test/chained-cylinder-spring.py	2010-10-20 11:31:49 +0000
+++ scripts/test/chained-cylinder-spring.py	2010-10-25 14:52:21 +0000
@@ -9,7 +9,7 @@
 poisson=5
 density=2.60e3 
 frictionAngle=radians(30)
-O.materials.append(CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,label='mat'))
+O.materials.append(CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,normalCohesion=1e13,shearCohesion=1e13,momentRotationLaw=True,label='mat'))
 O.dt=1e-4
 
 O.engines=[
@@ -20,8 +20,8 @@
 	]),
 	InteractionLoop(
 		[Ig2_ChainedCylinder_ChainedCylinder_ScGeom(),Ig2_Sphere_ChainedCylinder_CylScGeom()],
-		[Ip2_2xCohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True,normalCohesion=1e13,shearCohesion=1e13)],
-		[Law2_ScGeom_CohFrictPhys_CohesionMoment(momentRotationLaw=True,label='law')]
+		[Ip2_2xCohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True)],
+		[Law2_ScGeom_CohFrictPhys_CohesionMoment(label='law')]
 	),
 	## Apply gravity
 	GravityEngine(gravity=[0,-9.81,0]),