yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-03-09 21:57:46 +0100 (Mon, 09 Mar 2009)
New Revision: 1715

Modified:
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
   trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
   trunk/pkg/dem/SConscript
Log:
1. Added SpheresContactGeometry::updateShearForce, will be used (not activated though yet) by ElasticContactLaw and other.


Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2009-03-08 23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2009-03-09 20:57:46 UTC (rev 1715)
@@ -9,6 +9,62 @@
 // At least one virtual method must be in the .cpp file (!!!)
 SpheresContactGeometry::~SpheresContactGeometry(){};
 
+void SpheresContactGeometry::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, bool isDynamic1, bool isDynamic2, Real dt, bool avoidGranularRatcheting){
+
+	Vector3r axis;
+	Real angle;
+
+	// approximated rotations
+		axis = prevNormal.Cross(normal); 
+		shearForce -= shearForce.Cross(axis);
+		//angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angularVelocity+de2->angularVelocity);
+		//FIXME: if one body is kinematic then assumed its rotation centre does not lies along the normal
+		//(i.e. virtual sphere, which replaces this kinematic body on contact, does not rotate)
+		Vector3r summaryAngularVelocity(Vector3r::ZERO);
+		if (isDynamic1) summaryAngularVelocity += rbp1->angularVelocity;
+		if (isDynamic2) summaryAngularVelocity += rbp2->angularVelocity;
+		angle = dt*0.5*normal.Dot(summaryAngularVelocity);
+		axis = angle*normal;
+		shearForce -= shearForce.Cross(axis);
+		
+	// exact rotations
+	#if 0
+		Quaternionr q;
+		axis					= prevNormal.Cross(normal);
+		angle					= acos(normal.Dot(prevNormal));
+		q.FromAngleAxis(angle,axis);
+		shearForce        = shearForce*q;
+		angle             = dt*0.5*normal.dot(rbp1->angularVelocity+rbp2->angularVelocity);
+		axis					= normal;
+		q.FromAngleAxis(angle,axis);
+		shearForce        = q*shearForce;
+	#endif
+
+	Vector3r& x = contactPoint;
+	Vector3r c1x, c2x;
+
+	if(avoidGranularRatcheting){
+		/* The following definition of c1x and c2x is to avoid "granular ratcheting" 
+		 *  (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN, 
+		 *  "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
+		 *  ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors) */
+		c1x = (isDynamic1) ?  radius1*normal : x - rbp1->zeroPoint;
+		c2x = (isDynamic2) ? -radius2*normal : x - rbp2->zeroPoint;
+	}
+	else {
+		c1x = (x - rbp1->se3.position);
+		c2x = (x - rbp2->se3.position);
+	}
+
+	Vector3r relativeVelocity = (rbp2->velocity+rbp2->angularVelocity.Cross(c2x)) - (rbp1->velocity+rbp1->angularVelocity.Cross(c1x));
+	Vector3r shearVelocity = relativeVelocity-normal.Dot(relativeVelocity)*normal;
+	Vector3r shearDisplacement = shearVelocity*dt;
+	shearForce -= ks*shearDisplacement;
+}
+
+
+
+
 Vector3r SpheresContactGeometry::relRotVector() const{
 	Quaternionr relOri12=ori1.Conjugate()*ori2;
 	Quaternionr oriDiff=initRelOri12.Conjugate()*relOri12;

Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-03-08 23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-03-09 20:57:46 UTC (rev 1715)
@@ -5,6 +5,7 @@
 
 #include<yade/core/InteractionGeometry.hpp>
 #include<yade/lib-base/yadeWm3.hpp>
+#include<yade/pkg-common/RigidBodyParameters.hpp>
 /*! Class representing geometry of two spheres in contact.
  *
  * exactRot code
@@ -98,6 +99,9 @@
 
 		SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),facetContactFace(0.),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();}
 		virtual ~SpheresContactGeometry();
+
+		void updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const RigidBodyParameters* rbp1, const RigidBodyParameters* rbp2, bool isDynamic1, bool isDynamic2, Real dt, bool avoidGranularRatcheting=true);
+
 	REGISTER_ATTRIBUTES(/* no attributes from base class */,
 			(normal)
 			(contactPoint)

Modified: trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp	2009-03-08 23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/Engine/StandAloneEngine/ElasticContactLaw.cpp	2009-03-09 20:57:46 UTC (rev 1715)
@@ -114,6 +114,10 @@
 			Real un=currentContactGeometry->penetrationDepth;
 			currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
 	
+	#if 0
+		// the core under #else, refactored
+		currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,isDynamic1,isDynamic2,dt);
+	#else
 			Vector3r axis;
 			Real angle;
 
@@ -163,6 +167,8 @@
 			Vector3r shearVelocity			= relativeVelocity-currentContactGeometry->normal.Dot(relativeVelocity)*currentContactGeometry->normal;
 			Vector3r shearDisplacement		= shearVelocity*dt;
 			shearForce 			       -= currentContactPhysics->ks*shearDisplacement;
+
+	#endif
 	
 	// PFC3d SlipModel, is using friction angle. CoulombCriterion
 			Real maxFs = currentContactPhysics->normalForce.SquaredLength() * std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
@@ -173,7 +179,9 @@
 			}
 	////////// PFC3d SlipModel
 	
-			Vector3r f				= currentContactPhysics->normalForce + shearForce;
+			Vector3r f=currentContactPhysics->normalForce + shearForce;
+			Vector3r c1x(currentContactGeometry->contactPoint-de1->se3.position),
+				c2x(currentContactGeometry->contactPoint-de2->se3.position);
 			#ifdef BEX_CONTAINER
 				ncb->bex.addForce (id1,-f);
 				ncb->bex.addForce (id2,+f);

Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript	2009-03-08 23:03:42 UTC (rev 1714)
+++ trunk/pkg/dem/SConscript	2009-03-09 20:57:46 UTC (rev 1715)
@@ -54,15 +54,15 @@
 
 	env.SharedLibrary('SpheresContactGeometry',
 		['DataClass/InteractionGeometry/SpheresContactGeometry.cpp'],
-		LIBS=env['LIBS']+['yade-serialization','yade-base']),
+		LIBS=env['LIBS']+['RigidBodyParameters']),
 
 	env.SharedLibrary('SDECLinkGeometry',
-		['DataClass/InteractionGeometry/SDECLinkGeometry.cpp'],
-		LIBS=env['LIBS']+['yade-serialization','yade-base']),
+		['DataClass/InteractionGeometry/SDECLinkGeometry.cpp']
+	),
 
 	env.SharedLibrary('InteractionOfMyTetrahedron',
 		['DataClass/InteractionGeometry/InteractionOfMyTetrahedron.cpp'],
-		LIBS=env['LIBS']+['yade-serialization', 'yade-base']),
+	),
 
 	env.SharedLibrary('ElasticContactInteraction',
 		['DataClass/InteractionPhysics/ElasticContactInteraction.cpp'],