yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2008-09-17 23:48:17 +0200 (Wed, 17 Sep 2008)
New Revision: 1516

Modified:
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
   trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
   trunk/scripts/exact-rot.py
Log:
1. The exactRot code is reasonably verified now to be functional, including rolling correction (SpheresContactGeometry::relocateContactPoints) and will used in brefcom soon. 
2. Updated the "testing" script scripts/exact-rot.py, as a showcase as well.


Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2008-09-15 20:27:49 UTC (rev 1515)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.cpp	2008-09-17 21:48:17 UTC (rev 1516)
@@ -5,28 +5,31 @@
 #include "SpheresContactGeometry.hpp"
 YADE_PLUGIN("SpheresContactGeometry");
 
-/* Set contact points on both spheres such that their projection is the one given.
+/* Set contact points on both spheres such that their projection is the one given
+ * (should be on the plane passing through origin and oriented with normal; not checked!)
  */
 void SpheresContactGeometry::setTgPlanePts(Vector3r p1new, Vector3r p2new){
 	cp1rel=ori1.Conjugate()*rollPlanePtToSphere(p1new,d1,normal);
-	cp2rel=ori2.Conjugate()*rollPlanePtToSphere(-p2new,d2,-normal);
+	cp2rel=ori2.Conjugate()*rollPlanePtToSphere(p2new,d2,-normal);
 }
 
 
-/* Move contact point on both spheres in such way that their relative position is the same;
+/* Move contact point on both spheres in such way that their relative position (displacementT) is the same;
  * this should be done regularly to ensure that the angle doesn't go over π, since then quaternion would
  * flip axis and the point would project on other side of the tangent plane piece.
  */
-void SpheresContactGeometry::relocateContactPoint(){
+void SpheresContactGeometry::relocateContactPoints(){
 	Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
 	Vector3r midPt=.5*(p1+p2);
-	/* the factor 1.2 is arbitrary; it should be smaller than pi/2 and bigger than some reasonably small value to avoid frequent relocation */
-	if(midPt.Length()<1.2*min(d1,d2)) return;
-	setTgPlanePts(p1-midPt,p2-midPt);
+	if((p1.SquaredLength()>4*d1 || p2.SquaredLength()>4*d2) && midPt.SquaredLength()>.5*min(d1,d2)){
+		//cerr<<"RELOCATION with displacementT="<<displacementT(); // should be the same before and after relocation
+		setTgPlanePts(p1-midPt,p2-midPt);
+		//cerr<<" → "<<displacementT()<<endl;
+	}
 }
 
 /*! Perform slip of the projected contact points so that their distance becomes equal (or remains smaller) than the given one.
- *
+ * TODO: not yet tested
  */
 void SpheresContactGeometry::slipToDistIfNeeded(Real dist){
 	Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
@@ -49,8 +52,8 @@
  * @returns The projected point coordinates (with origin at the contact point).
  */
 Vector3r SpheresContactGeometry::unrollSpherePtToPlane(const Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& planeNormal){
-	Quaternionr fromNormalToPt; fromNormalToPt.Align(planeNormal,fromXtoPtOri*Vector3r::UNIT_X);
-	Vector3r axis; Real angle; fromNormalToPt.ToAxisAngle(axis,angle);
+	Quaternionr normal2pt; normal2pt.Align(planeNormal,fromXtoPtOri*Vector3r::UNIT_X);
+	Vector3r axis; Real angle; normal2pt.ToAxisAngle(axis,angle);
 	return (angle*radius) /* length */ *(axis.Cross(planeNormal)) /* direction: both are unit vectors */;
 }
 
@@ -68,7 +71,9 @@
 Quaternionr SpheresContactGeometry::rollPlanePtToSphere(const Vector3r& planePt, const Real& radius, const Vector3r& planeNormal){
 		Vector3r axis=planeNormal.Cross(planePt); axis.Normalize();
 		Real angle=planePt.Length()/radius;
-		return Quaternionr(axis,angle);
+		Quaternionr normal2pt(axis,angle);
+		Quaternionr ret; ret.Align(Vector3r::UNIT_X,normal2pt*planeNormal);
+		return ret;
 }
 
 

Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2008-09-15 20:27:49 UTC (rev 1515)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2008-09-17 21:48:17 UTC (rev 1516)
@@ -20,7 +20,9 @@
  * (a) plastic slip, when we want to limit their maximum distance (in the projected plane)
  * (b) rolling, where those points must be relocated to not flip over the π angle from the current contact point
  *
+ * TODO: add torsion code, that will (if a flag is on) compute torsion angle on the contact axis.
  *
+ *
  */
 class SpheresContactGeometry: public InteractionGeometry{
 	public :
@@ -39,28 +41,29 @@
 		 * 		taken instead)
 		 */
 		Quaternionr cp1rel, cp2rel;
-		// interaction "radii" and total length; this is _really_ contant throughout the interaction life
+		// interaction "radii" and total length; this is _really_ constant throughout the interaction life
 		// d1 is really distance from the sphere1 center to the contact plane, it may not correspond to the sphere radius!
 		// therefore, d1+d2=d0 (distance at which the contact was created)
 		Real d1, d2, d0;
 
-		// auxiliary functions: "this" not needed
+		// auxiliary functions; the quaternion magic is all in here
 		static Vector3r unrollSpherePtToPlane(const Quaternionr& fromXtoPtOri, const Real& radius, const Vector3r& normal);
 		static Quaternionr rollPlanePtToSphere(const Vector3r& planePt, const Real& radius, const Vector3r& normal);
+
 		void setTgPlanePts(Vector3r p1new, Vector3r p2new);
 
 		Vector3r contPtInTgPlane1(){return unrollSpherePtToPlane(ori1*cp1rel,d1,normal);}
 		Vector3r contPtInTgPlane2(){return unrollSpherePtToPlane(ori2*cp2rel,d2,-normal);}
 		Vector3r contPt(){return contactPoint; /*pos1+(d1/d0)*(pos2-pos1)*/}
 
-		Vector3r epsT(){return (1/d0)*(contPtInTgPlane2()-contPtInTgPlane1());}
-		Real epsN(){return (pos2-pos1).Length()/d0;}
+		Real displacementN(){return (pos2-pos1).Length()-d0;}
+		Real epsN(){return displacementN()*(1./d0);}
+		Vector3r displacementT(bool relocate=true){ return contPtInTgPlane2()-contPtInTgPlane1();}
+		Vector3r epsT(){return displacementT()*(1./d0);}
 	
 		void slipToDistIfNeeded(Real slip);
-		void relocateContactPoint();
+		void relocateContactPoints();
 
-
-
 		SpheresContactGeometry():InteractionGeometry(),contactPoint(Vector3r::ZERO),radius1(0),radius2(0),exactRot(false){createIndex();}
 		virtual ~SpheresContactGeometry(){};
 	protected :

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp	2008-09-15 20:27:49 UTC (rev 1515)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingSphere2InteractingSphere4SpheresContactGeometry.cpp	2008-09-17 21:48:17 UTC (rev 1516)
@@ -13,6 +13,7 @@
 #include<yade/pkg-common/InteractingSphere.hpp>
 
 #include<yade/lib-base/yadeWm3Extra.hpp>
+#include<yade/core/Omega.hpp>
 
 
 InteractingSphere2InteractingSphere4SpheresContactGeometry::InteractingSphere2InteractingSphere4SpheresContactGeometry()
@@ -36,17 +37,13 @@
 	InteractingSphere* s2 = static_cast<InteractingSphere*>(cm2.get());
 
 	Vector3r normal = se32.position-se31.position;
-	Real penetrationDepth = pow(interactionDetectionFactor*(s1->radius+s2->radius), 2) - normal.SquaredLength();// Compute a wrong value just to check sign - faster than computing distances
+	Real penetrationDepth = pow(interactionDetectionFactor*(s1->radius+s2->radius),2) - normal.SquaredLength();// Compute a wrong value just to check sign - faster than computing distances
 	//Real penetrationDepth = s1->radius+s2->radius-normal.Normalize();
 	if (penetrationDepth>0 || c->isReal){
 		shared_ptr<SpheresContactGeometry> scm;
-		if (c->interactionGeometry){
-			// WARNING! 
-			// FIXME - this must be dynamic cast until the contaners are rewritten to support multiple interactions types
-			//         the problem is that scm can be either SDECLinkGeometry or SpheresContactGeometry and the only way CURRENTLY
-			//         to check this is by dynamic cast. This has to be fixed.
-			scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
-		} else scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
+		if (c->interactionGeometry) scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
+		else scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
+
 		penetrationDepth = s1->radius+s2->radius-normal.Normalize();
 		scm->contactPoint = se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
 		scm->normal = normal;
@@ -58,9 +55,7 @@
 			scm->exactRot=true;
 			scm->pos1=se31.position; scm->pos2=se32.position;
 			scm->ori1=se31.orientation; scm->ori2=se32.orientation;
-			//scm->ori1.Normalize(); scm->ori2.Normalize();
 			if(c->isNew){
-				//cerr<<"+++ Assigning constants to SpheresContactGeometry"<<endl;
 				// contact constants
 				scm->d0=(se32.position-se31.position).Length();
 				scm->d1=s1->radius-penetrationDepth; scm->d2=s2->radius-penetrationDepth;
@@ -68,12 +63,9 @@
 				scm->cp1rel.Align(Vector3r::UNIT_X,se31.orientation.Conjugate()*normal);
 				scm->cp2rel.Align(Vector3r::UNIT_X,se32.orientation.Conjugate()*(-normal));
 				scm->cp1rel.Normalize(); scm->cp2rel.Normalize();
-				cerr<<"+++ Relative orientations: "<<scm->cp1rel<<" | "<<scm->cp2rel<<endl;
-				//cerr<<"+++ "<<se31.orientation.Conjugate()<<" | "<<se31.orientation.Conjugate()*normal<<"|"<<scm->cp1rel<<endl;
-				//cerr<<"@@@ cp1rel="<<scm->cp1rel[0]<<";"<<scm->cp1rel[1]<<";"<<scm->cp1rel[2]<<";"<<scm->cp1rel[3]<<endl;
-				//cerr<<"@@@ ori1="<<scm->ori1[0]<<";"<<scm->ori1[1]<<";"<<scm->ori1[2]<<";"<<scm->ori1[3]<<endl;
-				//cerr<<"+++ (normalized) "<<scm->cp1rel<<" || product with "<<se31.orientation<<" is "<<scm->cp1rel*se31.orientation<<endl;
 			}
+			// FIXME: temporary, testing only!
+			if((Omega::instance().getCurrentIteration()%500)==0) scm->relocateContactPoints();
 		}
 		return true;
 	} else return false;

Modified: trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp	2008-09-15 20:27:49 UTC (rev 1515)
+++ trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp	2008-09-17 21:48:17 UTC (rev 1516)
@@ -75,9 +75,9 @@
 	}
 
 	if(sc->exactRot){
-		//GLUtils::GLDrawLine(sc->pos1,sc->pos2,Vector3r(1,1,1));
+		GLUtils::GLDrawLine(sc->pos1,sc->pos2,Vector3r(.5,.5,.5));
 		// sphere center to point on the sphere
-		GLUtils::GLDrawText("[1]",sc->pos1,Vector3r(1,1,1)); GLUtils::GLDrawText("[2]",sc->pos2,Vector3r(1,1,1));
+		//GLUtils::GLDrawText("[1]",sc->pos1,Vector3r(1,1,1)); GLUtils::GLDrawText("[2]",sc->pos2,Vector3r(1,1,1));
 		GLUtils::GLDrawLine(sc->pos1,sc->pos1+(sc->ori1*sc->cp1rel*Vector3r::UNIT_X),Vector3r(0,.5,1));
 		GLUtils::GLDrawLine(sc->pos2,sc->pos2+(sc->ori2*sc->cp2rel*Vector3r::UNIT_X),Vector3r(0,1,.5));
 		//cerr<<"=== cp1rel="<<sc->cp1rel[0]<<";"<<sc->cp1rel[1]<<";"<<sc->cp1rel[2]<<";"<<sc->cp1rel[3]<<endl;
@@ -89,8 +89,7 @@
 		GLUtils::GLDrawLine(sc->contPt(),sc->contPt()+ptTg1,Vector3r(0,.5,1));
 		GLUtils::GLDrawLine(sc->contPt(),sc->contPt()+ptTg2,Vector3r(0,1,.5));
 		// projected shear
-		GLUtils::GLDrawLine(sc->contPt()+ptTg1,sc->contPt()+ptTg2,Vector3r(.7,.7,.7));
-		// TODO: crosshair to show contact plane (?)
+		GLUtils::GLDrawLine(sc->contPt()+ptTg1,sc->contPt()+ptTg2,Vector3r(1,1,1));
 	}
 
 

Modified: trunk/scripts/exact-rot.py
===================================================================
--- trunk/scripts/exact-rot.py	2008-09-15 20:27:49 UTC (rev 1515)
+++ trunk/scripts/exact-rot.py	2008-09-17 21:48:17 UTC (rev 1516)
@@ -22,7 +22,7 @@
 	StandAloneEngine('ElasticContactLaw'),
 	DeusExMachina('GravityEngine',{'gravity':[0,0,-9.81]}),
 	DeusExMachina('RotationEngine',{'subscribedBodies':[1],'rotationAxis':[1,0,0],'angularVelocity':.01}),
-	DeusExMachina('RotationEngine',{'subscribedBodies':[0],'rotationAxis':[sqrt(2)/2.,sqrt(2)/2.,0],'angularVelocity':.02}),
+	DeusExMachina('RotationEngine',{'subscribedBodies':[0],'rotationAxis':[1,1,1],'angularVelocity':-.02}),
 	MetaEngine('PhysicalActionDamper',[
 		EngineUnit('CundallNonViscousForceDamping',{'damping':0.2}),
 		EngineUnit('CundallNonViscousMomentumDamping',{'damping':0.2})
@@ -36,11 +36,11 @@
 ]
 from yade import utils
 o.bodies.append(utils.sphere([0,0,0],1,dynamic=False,color=[1,0,0],young=30e9,poisson=.3,density=2400,wire=True))
-o.bodies.append(utils.sphere([0,0,2],1,color=[0,1,0],young=30e9,poisson=.3,density=2400,wire=True))
+o.bodies.append(utils.sphere([0,sqrt(2),sqrt(2)],1,color=[0,1,0],young=30e9,poisson=.3,density=2400,wire=True))
 o.miscParams=[Generic('GLDrawSphere',{'glutUse':True})]
 
-o.dt=.2*utils.PWaveTimeStep()
-o.save('/tmp/a.xml.bz2');
+o.dt=.8*utils.PWaveTimeStep()
+#o.save('/tmp/a.xml.bz2');
 #o.run(100000); o.wait(); print o.iter/o.realtime,'iterations/sec'
 from yade import qt
 renderer=qt.Renderer()