← Back to team overview

yade-dev team mailing list archive

  1. yade-dev team
  2. Mailing list archive
  3. Message #06570

[Branch ~yade-dev/yade/trunk] Rev 2608: - recommit 2608 with fixed reg test (making ratcheting correction negligeable with dist=.99999999...

  Thread Previous • Date Previous • • Thread Next 
------------------------------------------------------------
revno: 2608
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Fri 2010-12-10 20:16:43 +0100
message:
  - recommit 2608 with fixed reg test (making ratcheting correction negligeable with dist=.999999999999999999*size)
  - also fix the compile errors
modified:
  pkg/common/Cylinder.cpp
  pkg/dem/ScGeom.cpp
  pkg/dem/ScGeom.hpp
  py/tests/pbc.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/common/Cylinder.cpp'
--- pkg/common/Cylinder.cpp	2010-12-06 19:24:20 +0000
+++ pkg/common/Cylinder.cpp	2010-12-10 19:16:43 +0000
@@ -65,6 +65,12 @@
 		direction = segment/length;
 		dist = direction.dot(branch);
 		if ((dist<-interactionDetectionFactor*cylinder->radius) && isNew) return false;
+		if (cylinderSt->rank>0){//make sure there is no contact with the previous element in the chain, or consider it a duplicate and continue
+			const shared_ptr<Body> cylinderPrev = Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank-1],scene);
+			Vector3r branchP = sphereSt->pos-cylinderPrev->state->pos;
+			Vector3r segmentP = cylinderSt->pos - cylinderPrev->state->pos;
+			if (segmentP.dot(branchP)<segmentP.dot(segmentP)) return false;//will give false on existing interaction, as expected
+		}
 	} else {//handle the last node with length=0
 		segment = Vector3r(0,0,0);
 		branch = sphereSt->pos-cylinderSt->pos-shift2;

=== modified file 'pkg/dem/ScGeom.cpp'
--- pkg/dem/ScGeom.cpp	2010-12-06 19:24:20 +0000
+++ pkg/dem/ScGeom.cpp	2010-12-10 19:16:43 +0000
@@ -37,37 +37,30 @@
 }
 
 Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shift2, const Vector3r& shiftVel, bool avoidGranularRatcheting){
-	Vector3r& x = contactPoint;
-	Vector3r c1x, c2x;
 	if(avoidGranularRatcheting){
 		//FIXME : put the long comment on the wiki and keep only a small abstract and link here.
-		/* B.C. Comment : The following definition of c1x and c2x is to avoid "granular ratcheting". This phenomenon has been introduced to me by S. McNamara in a seminar held in Paris, 2004 (GDR MiDi). The concept is also mentionned in many McNamara, Hermann, Lüding, and co-workers papers (see online : "Discrete element methods for the micro-mechanical investigation of granular ratcheting", R. García-Rojo, S. McNamara, H. J. Herrmann, Proceedings ECCOMAS 2004, @ http://www.ica1.uni-stuttgart.de/publications/2004/GMH04/), where it refers to an accumulation of strain in cyclic loadings.
-
-	        Unfortunately, published papers tend to focus on the "good" ratcheting, i.e. irreversible deformations due to the intrinsic nature of frictional granular materials, while the talk of McNamara in Paris clearly mentioned a possible "bad" ratcheting, purely linked with the formulation of the contact laws in what he called "molecular dynamics" (i.e. Cundall's model, as opposed to "contact dynamics" from Moreau and Jean).
-
-		Giving a short explanation :
-		The bad ratcheting is best understood considering this small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
+		/* B.C. Comment :
+		Giving a short explanation of what we want to avoid :
+		Numerical ratcheting is best understood considering a small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
 		1. translation "dx" in the normal direction
 		2. rotation "a"
 		3. translation "-dx" (back to initial position)
 		4. rotation "-a" (back to initial orientation)
 
-		If the branch vector used to define the relative shear in rotation*branch is not constant (typically if it is defined from the vector center->contactPoint like in the "else" below), then the shear displacement at the end of this cycle is not null : rotations a and -a are multiplied by branches of different lengths.
-		It results in a finite contact force at the end of the cycle even though the positions and orientations are unchanged, in total contradiction with the elastic nature of the problem. It could also be seen as an inconsistent energy creation or loss. It is BAD! And given the fact that DEM simulations tend to generate oscillations around equilibrium (damped mass-spring), it can have a significant impact on the evolution of the packings, resulting for instance in slow creep in iterations under constant load.
-		The solution to avoid that is quite simple : use a constant branch vector, like here radius_i*normal.
+		This problem has been analyzed in details by McNamara and co-workers. One will find interesting discussions in e.g. DOI 10.1103/PhysRevE.77.031304, even though solution it suggests is not fully applied here (equations of motion are not incorporating alpha, in contradiction with what is suggested by McNamara).
 		 */
-		// FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
-		c1x =   radius1*normal;
-		c2x =  -radius2*normal;
+		// For sphere-facet contact this will give an erroneous value of relative velocity...
+		Real alpha = (radius1+radius2)/(radius1+radius2-penetrationDepth);
+		Vector3r relativeVelocity = (rbp2->vel-rbp1->vel)*alpha + rbp2->angVel.cross(-radius2*normal) - rbp1->angVel.cross(radius1*normal);
+		relativeVelocity+=alpha*shiftVel;
+		return relativeVelocity;
 	} else {
-		// FIXME: It is correct for sphere-sphere and sphere-facet contact
-		c1x = (x - rbp1->pos);
-		c2x = (x - rbp2->pos + shift2);
-	}
-	Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
-	//update relative velocity for interactions across periods
-	relativeVelocity+=shiftVel;
-	return relativeVelocity;
+		// It is correct for sphere-sphere and sphere-facet contact
+		Vector3r c1x = (contactPoint - rbp1->pos);
+		Vector3r c2x = (contactPoint - rbp2->pos + shift2);
+		Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
+		relativeVelocity+=shiftVel;
+		return relativeVelocity;}
 }
 
 Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
@@ -90,10 +83,11 @@
 	if (isNew) {
 		initRotations(rbp1,rbp2);
 	} else {
-#ifdef YADE_SCGEOM_DEBUG
- 		rbp1.ori.normalize(); rbp2.ori.normalize(); initialOrientation2.normalize(); initialOrientation1.normalize();
-#endif
+// #ifdef YADE_SCGEOM_DEBUG
+ 		rbp1.ori.normalize(); rbp2.ori.normalize(); //initialOrientation2.normalize(); initialOrientation1.normalize();
+// #endif
 		Quaternionr delta((rbp1.ori * (initialOrientation1.conjugate())) * (initialOrientation2 * (rbp2.ori.conjugate())));
+		delta.normalize();
 		if (creep) delta = delta * twistCreep;
 		AngleAxisr aa(delta); // axis of rotation - this is the Moment direction UNIT vector; // angle represents the power of resistant ELASTIC moment
 		//Eigen::AngleAxisr(q) returns nan's when q close to identity, next tline fixes the pb.
@@ -123,10 +117,10 @@
 {
 	initialOrientation1	= state1.ori;
 	initialOrientation2	= state2.ori;
-	initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),normal);
-	currentContactOrientation = initialContactOrientation;
-	orientationToContact1   = initialOrientation1.conjugate() * initialContactOrientation;
-	orientationToContact2	= initialOrientation2.conjugate() * initialContactOrientation;
+// 	initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),normal);
+// 	currentContactOrientation = initialContactOrientation;
+// 	orientationToContact1   = initialOrientation1.conjugate() * initialContactOrientation;
+// 	orientationToContact2	= initialOrientation2.conjugate() * initialContactOrientation;
 	twist=0;
 	bending=Vector3r::Zero();
 	twistCreep=Quaternionr(1.0,0.0,0.0,0.0);

=== modified file 'pkg/dem/ScGeom.hpp'
--- pkg/dem/ScGeom.hpp	2010-11-07 11:46:20 +0000
+++ pkg/dem/ScGeom.hpp	2010-12-10 19:16:43 +0000
@@ -22,7 +22,7 @@
 		Vector3r twist_axis;//rotation vector around normal
 		Vector3r orthonormal_axis;//rotation vector in contact plane
 	public:
-		// inherited from GenericSpheresContact: Vector3r& normal; 
+		// inherited from GenericSpheresContact: Vector3r& normal;
 		Real &radius1, &radius2;
 		virtual ~ScGeom();
 		inline ScGeom& operator= (const ScGeom& source){
@@ -31,7 +31,7 @@
 			radius1=source.radius1; radius2=source.radius2;
 			penetrationDepth=source.penetrationDepth; shearInc=source.shearInc;
 			return *this;}
-		
+
 		//!precompute values of shear increment and interaction rotation data. Update contact normal to the currentNormal value. Precondition : the value of normal is not updated outside (and before) this function.
 		void precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r& currentNormal, bool isNew, const Vector3r& shift2, bool avoidGranularRatcheting=true);
 
@@ -43,10 +43,10 @@
 		Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, const Vector3r& shift2, bool avoidGranularRatcheting=true);
 		// Implement another version of getIncidentVel which does not handle periodicity.
 		Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
-		
+
 		// convenience version to be called from python
-		Vector3r getIncidentVel_py(shared_ptr<Interaction> i, bool avoidGranularRatcheting); 
-	
+		Vector3r getIncidentVel_py(shared_ptr<Interaction> i, bool avoidGranularRatcheting);
+
 		YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<IGeom>` of two :yref:`bodies<Body>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
 		((Real,penetrationDepth,NaN,(Attr::noSave|Attr::readonly),"Penetration distance of spheres (positive if overlapping)"))
 		((Vector3r,shearInc,Vector3r::Zero(),(Attr::noSave|Attr::readonly),"Shear displacement increment in the last step"))
@@ -67,15 +67,15 @@
 
 		void precomputeRotations(const State& rbp1, const State& rbp2, bool isNew, bool creep=false);
 		void initRotations(const State& rbp1, const State& rbp2);
-  		
+
 		YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom6D,ScGeom,"Class representing :yref:`geometry<IGeom>` of two :yref:`bodies<Body>` in contact. The contact has 6 DOFs (normal, 2×shear, twist, 2xbending) and uses :yref:`ScGeom` incremental algorithm for updating shear.",
-		((Quaternionr,orientationToContact1,Quaternionr(1.0,0.0,0.0,0.0),,""))
-		((Quaternionr,orientationToContact2,Quaternionr(1.0,0.0,0.0,0.0),,""))
+// 		((Quaternionr,orientationToContact1,Quaternionr(1.0,0.0,0.0,0.0),,""))
+// 		((Quaternionr,orientationToContact2,Quaternionr(1.0,0.0,0.0,0.0),,""))
 		((Quaternionr,initialOrientation1,Quaternionr(1.0,0.0,0.0,0.0),,""))
 		((Quaternionr,initialOrientation2,Quaternionr(1.0,0.0,0.0,0.0),,""))
 		((Quaternionr,twistCreep,Quaternionr(1.0,0.0,0.0,0.0),,""))
-		((Quaternionr,currentContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
-		((Quaternionr,initialContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
+// 		((Quaternionr,currentContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
+// 		((Quaternionr,initialContactOrientation,Quaternionr(1.0,0.0,0.0,0.0),,""))
 		((Real,twist,0,(Attr::noSave | Attr::readonly),"Elastic twist angle of the contact."))
 		((Vector3r,bending,Vector3r::Zero(),(Attr::noSave | Attr::readonly),"Bending at contact as a vector defining axis of rotation and angle (angle=norm)."))
 		,

=== modified file 'py/tests/pbc.py'
--- py/tests/pbc.py	2010-10-26 13:41:30 +0000
+++ py/tests/pbc.py	2010-12-10 19:16:43 +0000
@@ -20,7 +20,7 @@
 		O.reset(); O.periodic=True;
 		O.cell.refSize=Vector3(2.5,2.5,3)
 		self.cellDist=Vector3i(0,0,10) # how many cells away we go
-		self.relDist=Vector3(0,.999,0) # rel position of the 2nd ball within the cell
+		self.relDist=Vector3(0,.999999999999999999,0) # rel position of the 2nd ball within the cell
 		self.initVel=Vector3(0,0,5)
 		O.bodies.append(utils.sphere((1,1,1),.5))
 		self.initPos=Vector3([O.bodies[0].state.pos[i]+self.relDist[i]+self.cellDist[i]*O.cell.refSize[i] for i in (0,1,2)])
  Thread Previous • Date Previous • • Thread Next