yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2344
committer: bchareyre <bchareyre@dt-rv020>
branch nick: trunk
timestamp: Mon 2010-07-12 10:50:35 +0200
message:
  - miniEigen : memory leak source suspected (FIXME added)
  - ScGeom : #define IGCACHE by default, and keep behaviour unchanged for older functions (still some cleaning to do in this class, after all Law2's have been adapted).
  - CohFrictLaw : Eigen is returning nan quaternions after trivial operations, it is workarounded in the law.
  - cylinder : a new shape, first step in experiments for  tracking interactions jumping over chained elements.
  - utils.py : a demo script for chained cylinders  
added:
  pkg/common/DataClass/Shape/Cylinder.cpp
  pkg/common/DataClass/Shape/Cylinder.hpp
  scripts/test/chained-cylinder-spring.py
modified:
  pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp
  pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp
  pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
  py/mathWrap/miniEigen.cpp
  py/utils.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

=== added file 'pkg/common/DataClass/Shape/Cylinder.cpp'
--- pkg/common/DataClass/Shape/Cylinder.cpp	1970-01-01 00:00:00 +0000
+++ pkg/common/DataClass/Shape/Cylinder.cpp	2010-07-12 08:50:35 +0000
@@ -0,0 +1,257 @@
+#include "Cylinder.hpp"
+#include<yade/pkg-common/Sphere.hpp>
+#include<yade/lib-opengl/OpenGLWrapper.hpp>
+#include<yade/pkg-common/Aabb.hpp>
+
+Cylinder::~Cylinder(){}
+ChainedCylinder::~ChainedCylinder(){}
+ChainedState::~ChainedState(){}
+CylScGeom::~CylScGeom(){}
+
+
+YADE_PLUGIN((Cylinder)(ChainedCylinder)(ChainedState)(CylScGeom)(Ig2_Sphere_ChainedCylinder_CylScGeom)(Ig2_ChainedCylinder_ChainedCylinder_ScGeom)(Gl1_Cylinder)(Bo1_Cylinder_Aabb)/*(Bo1_ChainedCylinder_Aabb)*/);
+YADE_REQUIRE_FEATURE(OPENGL)
+
+vector<vector<int> > ChainedState::chains;
+unsigned int ChainedState::currentChain=0;
+
+//!##################	IG FUNCTORS   #####################
+
+#ifdef YADE_DEVIRT_FUNCTORS
+bool Ig2_Sphere_ChainedCylinder_CylScGeom::go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){ throw runtime_error("Do not call Ig2_Sphere_ChainedCylinder_CylScGeom::go, use getStaticFunctorPtr and call that function instead."); }
+bool Ig2_Sphere_ChainedCylinder_CylScGeom::goStatic(InteractionGeometryFunctor* _self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
+	const Ig2_Sphere_ChainedCylinder_CylScGeom* self=static_cast<Ig2_Sphere_ChainedCylinder_CylScGeom*>(_self);
+	const Real& interactionDetectionFactor=self->interactionDetectionFactor;
+#else
+bool Ig2_Sphere_ChainedCylinder_CylScGeom::go(	const shared_ptr<Shape>& cm1,
+							const shared_ptr<Shape>& cm2,
+							const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
+							const shared_ptr<Interaction>& c)
+{
+#endif
+	const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
+
+	Sphere *s1=static_cast<Sphere*>(cm1.get()), *s2=static_cast<Sphere*>(cm2.get());
+	Vector3r normal=(se32.position+shift2)-se31.position;
+	Real penetrationDepthSq=pow(interactionDetectionFactor*(s1->radius+s2->radius),2) - normal.squaredNorm();
+	if (penetrationDepthSq>0 || c->isReal() || force){
+		shared_ptr<ScGeom> scm;
+		bool isNew = !c->interactionGeometry;
+		if(!isNew) scm=YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
+		else { scm=shared_ptr<ScGeom>(new ScGeom()); c->interactionGeometry=scm; }
+
+		Real norm=normal.norm(); normal/=norm; // normal is unit vector now
+		Real penetrationDepth=s1->radius+s2->radius-norm;
+		scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
+		if(isNew) scm->prevNormal=normal;
+		else scm->prevNormal=scm->normal;
+		scm->normal=normal;
+		scm->penetrationDepth=penetrationDepth;
+		scm->radius1=s1->radius;
+		scm->radius2=s2->radius;
+#ifdef IGCACHE
+		if (scene->isPeriodic) {
+			Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
+ 			scm->precompute(state1,state2,scene->dt,shiftVel,true);}
+ 		else scm->precompute(state1,state2,scene->dt,true);
+#endif
+		return true;
+	}
+	return false;
+}
+
+
+bool Ig2_Sphere_ChainedCylinder_CylScGeom::goReverse(	const shared_ptr<Shape>& cm1,
+								const shared_ptr<Shape>& cm2,
+								const State& state1,
+								const State& state2,
+								const Vector3r& shift2,
+								const bool& force,
+								const shared_ptr<Interaction>& c)
+{
+	return go(cm1,cm2,state2,state1,-shift2,force,c);
+}
+
+
+#ifdef YADE_DEVIRT_FUNCTORS
+bool Ig2_ChainedCylinder_ChainedCylinder_ScGeom::go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){ throw runtime_error("Do not call Ig2_Sphere_ChainedCylinder_CylScGeom::go, use getStaticFunctorPtr and call that function instead."); }
+bool Ig2_ChainedCylinder_ChainedCylinder_ScGeom::goStatic(InteractionGeometryFunctor* _self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c){
+	const Ig2_ChainedCylinder_ChainedCylinder_ScGeom* self=static_cast<Ig2_ChainedCylinder_ChainedCylinder_ScGeom*>(_self);
+	const Real& interactionDetectionFactor=self->interactionDetectionFactor;
+#else
+bool Ig2_ChainedCylinder_ChainedCylinder_ScGeom::go(	const shared_ptr<Shape>& cm1,
+							const shared_ptr<Shape>& cm2,
+							const State& state1, const State& state2, const Vector3r& shift2, const bool& force,
+							const shared_ptr<Interaction>& c)
+{
+#endif
+	const ChainedState *pChain1, *pChain2;
+	pChain1=YADE_CAST<const ChainedState*>(&state1);
+	pChain2=YADE_CAST<const ChainedState*>(&state2);
+	const ChainedState& chain1 = *pChain1;
+	const ChainedState& chain2 = *pChain2;
+	if (!pChain1 || !pChain2) {
+		cerr <<"cast failed8567"<<endl;
+	}
+	bool revert = (pChain2->rank-chain1.rank == -1);
+
+	if (pChain2->chainNumber!=pChain1->chainNumber) {cerr<<"PROBLEM0124"<<endl; return false;}
+	else if (chain2.rank-chain1.rank != 1 && pChain2->rank-pChain1->rank != -1) {/*cerr<<"PROBLEM0123"<<endl;*/ return false;}
+	const ChainedState* pChain = revert ? pChain2 : pChain1;
+
+	ChainedCylinder *s1=static_cast<ChainedCylinder*>(cm1.get()), *s2=static_cast<ChainedCylinder*>(cm2.get());
+	ChainedCylinder *s = revert ? s2 : s1;
+	shared_ptr<ScGeom> scm;
+	bool isNew = !c->interactionGeometry;
+	if(!isNew) scm=YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
+	else { scm=shared_ptr<ScGeom>(new ScGeom()); c->interactionGeometry=scm; }
+	Real length=abs((pChain->ori*Vector3r::UnitZ()).dot(chain2.pos-chain1.pos));
+ 	Vector3r segment =pChain->ori*Vector3r::UnitZ()*length;
+	
+	if (revert) segment = -segment; 
+	if(isNew) {scm->normal=scm->prevNormal=segment/length;s->initLength=length;}
+	else {scm->prevNormal=scm->normal; scm->normal=segment/length;}
+
+	if (!revert) {
+		scm->radius1=s->initLength;
+		scm->radius2=0;}
+	else {
+		scm->radius1=0;
+		scm->radius2=s->initLength;}
+		
+	s->length=length;
+	scm->penetrationDepth=s->initLength-length;
+	scm->contactPoint=pChain->pos+pChain->ori*Vector3r::UnitZ()*length;
+#ifdef IGCACHE
+	if (scene->isPeriodic) {
+		Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
+		scm->precompute(state1,state2,scene->dt,shiftVel,true);}
+	else {
+		scm->precompute(state1,state2,scene->dt,true);
+	}
+#else
+	cerr<<"this is not supposed to work without #define IGCACHE"<<endl;
+#endif
+	//Set values that will be considered in Ip2 functor, geometry (precomputed) is really defined with values above
+	scm->radius1=s->initLength*0.5;
+	scm->radius2=s->initLength*0.5;
+	return true;
+}
+
+
+bool Ig2_ChainedCylinder_ChainedCylinder_ScGeom::goReverse(	const shared_ptr<Shape>& cm1,
+								const shared_ptr<Shape>& cm2,
+								const State& state1,
+								const State& state2,
+								const Vector3r& shift2,
+								const bool& force,
+								const shared_ptr<Interaction>& c)
+{
+	return go(cm2,cm1,state2,state1,-shift2,force,c);
+}
+
+//!##################	RENDERING   #####################
+
+bool Gl1_Cylinder::wire;
+bool Gl1_Cylinder::glutNormalize;
+int  Gl1_Cylinder::glutSlices;
+int  Gl1_Cylinder::glutStacks;
+int Gl1_Cylinder::glCylinderList=-1;
+
+void Gl1_Cylinder::go(const shared_ptr<Shape>& cm, const shared_ptr<State>& ,bool wire2, const GLViewInfo&)
+{
+	Real r=(static_cast<Cylinder*>(cm.get()))->radius;
+	Real length=(static_cast<Cylinder*>(cm.get()))->length;
+	//glMaterialv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, Vector3f(cm->color[0],cm->color[1],cm->color[2]));
+	glColor3v(cm->color);
+	if(glutNormalize)	glPushAttrib(GL_NORMALIZE); // as per http://lists.apple.com/archives/Mac-opengl/2002/Jul/msg00085.html
+// 	glPushMatrix();
+	if (wire || wire2) drawCylinder(true, r,length);
+	else drawCylinder(false, r,length);
+	if(glutNormalize) glPopAttrib();
+// 	glPopMatrix();
+	return;
+}
+
+
+void Gl1_Cylinder::drawCylinder(bool wire, Real radius, Real length) const
+{
+//    GLERROR;
+/*	if (glCylinderList<0) {
+		glCylinderList = glGenLists(1);
+		glNewList(glCylinderList,GL_COMPILE);*/
+   glPushMatrix();
+   GLUquadricObj *quadObj = gluNewQuadric();
+   gluQuadricDrawStyle(quadObj, (GLenum) (wire ? GLU_SILHOUETTE : GLU_FILL));
+   gluQuadricNormals(quadObj, (GLenum) GLU_SMOOTH);
+   gluQuadricOrientation(quadObj, (GLenum) GLU_OUTSIDE);
+//     glTranslatef(0.0,0.0,-length*0.5);
+   //scaling needs to adapt spheres or they will be elipsoids. They actually seem to disappear when commented glList code is uncommented, the cylinders are displayed correclty.
+//    glScalef(1,length,1);
+   gluCylinder(quadObj, radius, radius, length, glutSlices,glutStacks);
+   gluQuadricOrientation(quadObj, (GLenum) GLU_INSIDE);
+   glutSolidSphere(radius,glutSlices,glutStacks);
+   glTranslatef(0.0,0.0,length);
+//    glRotatef(180,0.0,1.0,0.0);
+   glutSolidSphere(radius,glutSlices,glutStacks);
+//    gluDisk(quadObj,0.0,radius,glutSlices,_loops);
+   gluDeleteQuadric(quadObj);
+   glPopMatrix();
+//    GLERROR;
+	
+// 	glEndList();}
+// 	glCallList(glCylinderList);
+
+}
+
+//!##################	BOUNDS FUNCTOR   #####################
+
+void Bo1_Cylinder_Aabb::go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r& se3, const Body* b){
+	Cylinder* cylinder = static_cast<Cylinder*>(cm.get());
+	Aabb* aabb = static_cast<Aabb*>(bv.get());
+	if(!scene->isPeriodic){
+		const Vector3r& O = se3.position;
+		Vector3r O2 = se3.position+se3.orientation*cylinder->segment;
+		aabb->min=aabb->max=O;
+		for (int k=0;k<3;k++){
+			aabb->min[k]=min(aabb->min[k],min(O[k],O2[k])-cylinder->radius);
+			aabb->max[k]=max(aabb->max[k],max(O[k],O2[k])+cylinder->radius);
+		}
+		return;
+	}
+	// adjust box size along axes so that cylinder doesn't stick out of the box even if sheared (i.e. parallelepiped)
+// 	if(scene->cell->hasShear()) {
+// 		Vector3r refHalfSize(minkSize);
+// 		const Vector3r& cos=scene->cell->getCos();
+// 		for(int i=0; i<3; i++){
+// 			//cerr<<"cos["<<i<<"]"<<cos[i]<<" ";
+// 			int i1=(i+1)%3,i2=(i+2)%3;
+// 			minkSize[i1]+=.5*refHalfSize[i1]*(1/cos[i]-1);
+// 			minkSize[i2]+=.5*refHalfSize[i2]*(1/cos[i]-1);
+// 		}
+// 	}
+// 	//cerr<<" || "<<halfSize<<endl;
+// 	aabb->min = scene->cell->unshearPt(se3.position)-minkSize;
+// 	aabb->max = scene->cell->unshearPt(se3.position)+minkSize;
+}
+/*
+void Bo1_ChainedCylinder_Aabb::go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r& se3, const Body* b){
+	Cylinder* cylinder = static_cast<Cylinder*>(cm.get());
+	Aabb* aabb = static_cast<Aabb*>(bv.get());
+	if(!scene->isPeriodic){
+		const Vector3r& O = se3.position-se3.positioncylinder->segment*0.5;
+		//In this case, segment will be updated by Ig2_ChainedCylinder_ChainedCylinder_ScGeom, no need to rotate it
+		Vector3r O2 = se3.position+cylinder->segment*0.5;
+// 		cerr << O<<" "<<O2<<endl;
+		aabb->min=aabb->max=O;
+		for (int k=0;k<3;k++){
+			aabb->min[k]=min(aabb->min[k],min(O[k],O2[k])-cylinder->radius);
+			aabb->max[k]=max(aabb->max[k],max(O[k],O2[k])+cylinder->radius);
+		}
+		return;
+	}
+}*/
+
+
+
+

=== added file 'pkg/common/DataClass/Shape/Cylinder.hpp'
--- pkg/common/DataClass/Shape/Cylinder.hpp	1970-01-01 00:00:00 +0000
+++ pkg/common/DataClass/Shape/Cylinder.hpp	2010-07-12 08:50:35 +0000
@@ -0,0 +1,167 @@
+#pragma once
+#include<yade/core/Shape.hpp>
+#include<yade/core/State.hpp>
+#include<yade/core/Body.hpp>
+#include<yade/pkg-dem/ScGeom.hpp>
+#include<yade/pkg-common/InteractionGeometryFunctor.hpp>
+#include<yade/core/Scene.hpp>
+#include<yade/pkg-common/GLDrawFunctors.hpp>
+#include<yade/pkg-common/BoundFunctor.hpp>
+
+
+class Cylinder: public Shape{
+	public:
+		Cylinder(Real _radius, Real _length): radius(_radius), length(_length), segment(Vector3r(0,0,1)*_length){ createIndex(); }
+		virtual ~Cylinder ();
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(Cylinder,Shape,"Geometry of a cylinder, as Minkowsky sum of line and sphere.",
+		((Real,radius,NaN,"Radius [m]"))
+		((Real,length,NaN,"Length [m]"))
+		((Vector3r,segment,Vector3r::Zero(),"Length vector")),
+		createIndex(); /*ctor*/
+	);
+	REGISTER_CLASS_INDEX(Cylinder,Shape);
+};
+
+class ChainedCylinder: public Cylinder{
+	public:
+		ChainedCylinder(Real _radius, Real _length): Cylinder(_radius,_length){}
+		virtual ~ChainedCylinder ();
+		//Keep pointers or copies of connected states?
+// 		ChainedState st1, st2;
+		
+
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(ChainedCylinder,Cylinder,"Geometry of a deformable chained cylinder, using geometry :yref:`MinkCylinder`.",
+  		((Real,initLength,0,"tensile-free length, used as reference for tensile strain"))
+		,createIndex();/*ctor*/
+// 		state=shared_ptr<ChainedState>(new ChainedState);
+
+	);
+	REGISTER_CLASS_INDEX(ChainedCylinder,Cylinder);
+};
+
+class CylScGeom: public ScGeom{
+	public:
+		virtual ~CylScGeom ();
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CylScGeom,ScGeom,"Geometry of a cylinder-sphere contact.",
+		((bool,onNode,false,"contact on node")),
+		createIndex(); /*ctor*/
+	);
+	REGISTER_CLASS_INDEX(CylScGeom,ScGeom);
+};
+
+
+class ChainedState: public State{
+	public:
+		static vector<vector<int> > chains;
+		static unsigned int currentChain;
+		vector<int> barContacts;
+		vector<int> nodeContacts;
+
+		virtual ~ChainedState ();
+		void addToChain(int bodyId) {
+			if (chains.size()<=currentChain) chains.resize(currentChain+1);
+			chainNumber=currentChain;
+ 			rank=chains[currentChain].size();
+ 			chains[currentChain].push_back(rank);}
+
+	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ChainedState,State,"State of a chained bodies, containing information on connectivity in order to track contacts jumping over contiguous elements. Chains are 1D lists from which id of chained bodies are retrieved via :yref:rank<ChainedState::rank>` and :yref:chainNumber<ChainedState::chainNumber>`.",
+ 		((int,rank,0,"rank in the chain"))
+ 		((int,chainNumber,0,"chain id"))
+		,
+		/* additional initializers */
+/*			((pos,se3.position))
+			((ori,se3.orientation)),*/
+		/* ctor */,
+		/*py*/,
+// 		.def_readwrite("chains",&ChainedState::chains,"documentation")
+		.def_readwrite("currentChain",&ChainedState::currentChain,"Current active chain (where newly created chained bodies will be appended).")
+		.def("addToChain",&ChainedState::addToChain,(python::arg("bodyId")),"Add body to current active chain")
+	);
+};
+
+
+
+class Ig2_Sphere_ChainedCylinder_CylScGeom: public InteractionGeometryFunctor{
+	public:
+		virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+		virtual bool goReverse(	const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+	#ifdef YADE_DEVIRT_FUNCTORS
+		void* getStaticFuncPtr(){ return (void*)&Ig2_Sphere_ChainedCylinder_CylScGeom::goStatic; }
+		static bool goStatic(InteractionGeometryFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+	#endif
+	YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_ChainedCylinder_CylScGeom,InteractionGeometryFunctor,"Create/update a :yref:`ScGeom` instance representing intersection of two :yref:`Spheres<Sphere>`.",
+		((Real,interactionDetectionFactor,1,"Enlarge both radii by this factor (if >1), to permit creation of distant interactions."))
+	);
+	FUNCTOR2D(Sphere,ChainedCylinder);
+	DEFINE_FUNCTOR_ORDER_2D(Sphere,ChainedCylinder);
+};
+
+class Ig2_ChainedCylinder_ChainedCylinder_ScGeom: public InteractionGeometryFunctor{
+	public:
+		virtual bool go(const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+		virtual bool goReverse(	const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+	#ifdef YADE_DEVIRT_FUNCTORS
+		void* getStaticFuncPtr(){ return (void*)&Ig2_Sphere_ChainedCylinder_CylScGeom::goStatic; }
+		static bool goStatic(InteractionGeometryFunctor* self, const shared_ptr<Shape>& cm1, const shared_ptr<Shape>& cm2, const State& state1, const State& se32, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& c);
+	#endif
+	YADE_CLASS_BASE_DOC_ATTRS(Ig2_ChainedCylinder_ChainedCylinder_ScGeom,InteractionGeometryFunctor,"Create/update a :yref:`ScGeom` instance representing connexion between :yref:`chained cylinders<ChainedCylinder>`.",
+		((Real,interactionDetectionFactor,1,"Enlarge both radii by this factor (if >1), to permit creation of distant interactions."))
+	);
+	FUNCTOR2D(ChainedCylinder,ChainedCylinder);
+	// needed for the dispatcher, even if it is symmetric
+	DEFINE_FUNCTOR_ORDER_2D(ChainedCylinder,ChainedCylinder);
+};
+
+
+
+
+class Gl1_Cylinder : public GlShapeFunctor{
+	private:
+		static int glCylinderList;
+		void subdivideTriangle(Vector3r& v1,Vector3r& v2,Vector3r& v3, int depth);
+		void drawCylinder(bool wire, Real radius, Real length) const;
+		void initGlLists(void);
+	public:
+		virtual void go(const shared_ptr<Shape>&, const shared_ptr<State>&,bool,const GLViewInfo&);
+	YADE_CLASS_BASE_DOC_STATICATTRS(Gl1_Cylinder,GlShapeFunctor,"Renders :yref:`Cylinder` object",
+		((bool,wire,false,"Only show wireframe (controlled by ``glutSlices`` and ``glutStacks``."))
+		((bool,glutNormalize,true,"Fix normals for non-wire rendering"))
+		((int,glutSlices,8,"Number of sphere slices."))
+		((int,glutStacks,4,"Number of sphere stacks."))
+	);
+	RENDERS(Cylinder);
+};
+
+class Bo1_Cylinder_Aabb : public BoundFunctor
+{
+	public :
+		void go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r&, const Body*);
+	FUNCTOR2D(Cylinder,Aabb);
+	YADE_CLASS_BASE_DOC_ATTRS(Bo1_Cylinder_Aabb,BoundFunctor,"Functor creating :yref:`Aabb` from :yref:`Cylinder`.",
+		((Real,aabbEnlargeFactor,((void)"deactivated",-1),"Relative enlargement of the bounding box; deactivated if negative.\n\n.. note::\n\tThis attribute is used to create distant interaction, but is only meaningful with an :yref:`InteractionGeometryFunctor` which will not simply discard such interactions: :yref:`Ig2_Cylinder_Cylinder_Dem3DofGeom::distFactor` / :yref:`Ig2_Cylinder_Cylinder_ScGeom::interactionDetectionFactor` should have the same value as :yref:`aabbEnlargeFactor<Bo1_Cylinder_Aabb::aabbEnlargeFactor>`."))
+	);
+};
+
+// Keep this : Cylinders and ChainedCylinders will have different centers maybe.
+// class Bo1_ChainedCylinder_Aabb : public Bo1_Cylinder_Aabb
+// {
+// 	public :
+// 		void go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r&, const Body*);
+// 	FUNCTOR2D(ChainedCylinder,Aabb);
+// 	YADE_CLASS_BASE_DOC_ATTRS(Bo1_ChainedCylinder_Aabb,Bo1_Cylinder_Aabb,"Functor creating :yref:`Aabb` from :yref:`Cylinder`.",
+// 		((Real,aabbEnlargeFactor,((void)"deactivated",-1),"Relative enlargement of the bounding box; deactivated if negative.\n\n.. note::\n\tThis attribute is used to create distant interaction, but is only meaningful with an :yref:`InteractionGeometryFunctor` which will not simply discard such interactions: :yref:`Ig2_Cylinder_Cylinder_Dem3DofGeom::distFactor` / :yref:`Ig2_Cylinder_Cylinder_ScGeom::interactionDetectionFactor` should have the same value as :yref:`aabbEnlargeFactor<Bo1_Cylinder_Aabb::aabbEnlargeFactor>`."))
+// 	);
+// };
+
+
+
+
+REGISTER_SERIALIZABLE(Bo1_Cylinder_Aabb);
+// REGISTER_SERIALIZABLE(Bo1_ChainedCylinder_Aabb);
+REGISTER_SERIALIZABLE(Gl1_Cylinder);
+REGISTER_SERIALIZABLE(Cylinder);
+REGISTER_SERIALIZABLE(ChainedCylinder);
+REGISTER_SERIALIZABLE(ChainedState);
+REGISTER_SERIALIZABLE(CylScGeom);
+REGISTER_SERIALIZABLE(Ig2_Sphere_ChainedCylinder_CylScGeom);
+REGISTER_SERIALIZABLE(Ig2_ChainedCylinder_ChainedCylinder_ScGeom);

=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2010-07-08 18:37:01 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp	2010-07-12 08:50:35 +0000
@@ -126,10 +126,10 @@
 
 //Removing this function will need to adapt all laws
 Vector3r ScGeom::rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
-#ifdef IGCACHE
-	rotate(shearForce);
-	return shearIncrement;
-#else
+// #ifdef IGCACHE
+// 	rotate(shearForce);
+// 	return shearIncrement;
+// #else
 	rotate(shearForce,prevNormal,rbp1,rbp2,dt);
 	Vector3r& x = contactPoint;
 	Vector3r c1x, c2x;
@@ -150,13 +150,13 @@
 	relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
 	Vector3r shearDisplacement = relativeVelocity*dt;
 	return shearDisplacement;
-#endif //IGCACHE
+// #endif //IGCACHE
 }
 
 Vector3r& ScGeom::rotate(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt){
-#ifdef IGCACHE
-	return rotate(shearForce);
-#else
+// #ifdef IGCACHE
+// 	return rotate(shearForce);
+// #else
 	Vector3r axis;
 	Real angle;
 	// approximated rotations
@@ -178,7 +178,7 @@
 // 		shearForce        = q*shearForce;
 
 	return shearForce;
-#endif
+// #endif
 }
 
 Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){

=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-07-08 18:37:01 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp	2010-07-12 08:50:35 +0000
@@ -14,7 +14,7 @@
  */
 
 #define SCG_SHEAR
-// #define IGCACHE
+#define IGCACHE
 
 class ScGeom: public GenericSpheresContact {
 	public:

=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp	2010-06-08 22:25:00 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp	2010-07-12 08:50:35 +0000
@@ -74,6 +74,7 @@
 	if (un < 0 && (currentContactPhysics->normalForce.squaredNorm() > pow(currentContactPhysics->normalAdhesion,2)
 	               || currentContactPhysics->normalAdhesion==0)) {
 		// BREAK due to tension
+		cerr <<"requesterase";
 		ncb->interactions->requestErase(contact->getId1(),contact->getId2());
 		// contact->interactionPhysics was reset now; currentContactPhysics still hold the object, but is not associated with the interaction anymore
 // 			currentContactPhysics->cohesionBroken = true;
@@ -85,7 +86,12 @@
 		///////////////////////// CREEP START ///////////
 		if (shear_creep) shearForce -= currentContactPhysics->ks*(shearForce*dt/creep_viscosity);
 		///////////////////////// CREEP END ////////////
+#ifdef IGCACHE
+		Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
+		const Vector3r& dus = currentContactGeometry->shearIncrement;
+#else
 		Vector3r dus = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt);
+#endif
 		//Linear elasticity giving "trial" shear force
 		shearForce -= currentContactPhysics->ks*dus;
 
@@ -98,13 +104,13 @@
 			if (currentContactPhysics->fragile && !currentContactPhysics->cohesionBroken) {
 				currentContactPhysics->SetBreakingState();
 				maxFs = max((Real) 0, Fn*currentContactPhysics->tangensOfFrictionAngle);
+				
 			}
 			maxFs = maxFs / Fs;
 			if (maxFs>1) cerr << "maxFs>1!!" << endl;
 			shearForce *= maxFs;
 			if (Fn<0)  currentContactPhysics->normalForce = Vector3r::Zero();//Vector3r::Zero()
 		}
-
 		applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
 
 		/// Moment law        ///
@@ -116,16 +122,16 @@
 			// Quaternionr align(axis,angle);
 			// currentContactPhysics->currentContactOrientation =  align * currentContactPhysics->currentContactOrientation;
 			//}
-
-			Quaternionr delta(b1->state->ori * currentContactPhysics->initialOrientation1.conjugate() *
-			                  currentContactPhysics->initialOrientation2 * b2->state->ori.conjugate());
+			Quaternionr delta((b1->state->ori * (currentContactPhysics->initialOrientation1.conjugate())) *
+			                  (currentContactPhysics->initialOrientation2 * (b2->state->ori.conjugate())));
 			if (twist_creep) {
 				delta = delta * currentContactPhysics->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.
+			if (isnan(aa.angle())) aa.angle()=0;
 			if (aa.angle() > Mathr::PI) aa.angle() -= Mathr::TWO_PI;   // angle is between 0 and 2*pi, but should be between -pi and pi
-
 			Real angle_twist(aa.angle() * aa.axis().dot(currentContactGeometry->normal));
 			Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
 

=== modified file 'py/mathWrap/miniEigen.cpp'
--- py/mathWrap/miniEigen.cpp	2010-06-08 22:25:00 +0000
+++ py/mathWrap/miniEigen.cpp	2010-07-12 08:50:35 +0000
@@ -99,6 +99,7 @@
 static Matrix3r* Matrix3r_fromElements(Real m00, Real m01, Real m02, Real m10, Real m11, Real m12, Real m20, Real m21, Real m22){ Matrix3r* m(new Matrix3r); (*m)<<m00,m01,m02,m10,m11,m12,m20,m21,m22; return m; }
 static Quaternionr Quaternionr_setFromTwoVectors(Quaternionr& q, const Vector3r& u, const Vector3r& v){ return q.setFromTwoVectors(u,v); }
 static Vector3r Quaternionr_Rotate(Quaternionr& q, const Vector3r& u){ return q*u; }
+//FIXME : risk of memory leak with these "new"?
 static Quaternionr* Quaternionr_fromAxisAngle(const Vector3r& axis, const Real angle){ return new Quaternionr(AngleAxisr(angle,axis)); }
 static Quaternionr* Quaternionr_fromAngleAxis(const Real angle, const Vector3r& axis){ return new Quaternionr(AngleAxisr(angle,axis)); }
 static bp::tuple Quaternionr_toAxisAngle(const Quaternionr& self){ AngleAxisr aa(self); return bp::make_tuple(aa.axis(),aa.angle());}

=== modified file 'py/utils.py'
--- py/utils.py	2010-06-29 21:02:29 +0000
+++ py/utils.py	2010-07-12 08:50:35 +0000
@@ -208,6 +208,36 @@
 	b.mask=mask
 	return b
 
+
+def chCylinder(begin=Vector3(0,0,0),end=Vector3(1.,0.,0.),radius=0.2,dynamic=True,wire=False,color=None,highlight=False,material=-1,mask=1):
+	"""Create and chain a MinkCylinder with given parameters. This shape is the Minkowsky sum of line and sphere.
+
+	:Parameters:
+		`radius`: Real
+			radius of sphere in the Minkowsky sum.
+		`begin`: Vector3
+			first point positioning the line in the Minkowsky sum
+		`last`: Vector3
+			last point positioning the line in the Minkowsky sum
+
+	In order to build a correct chain, last point of element of rank N must correspond to first point of element of rank N+1 in the same chain (with some tolerance, since bounding boxes will be used to create connections."""
+	segment=end-begin
+	b=Body()
+	b.shape=ChainedCylinder(radius=radius,length=segment.norm(),color=color if color else randomColor(),wire=wire,highlight=highlight)
+	b.shape.segment=b.shape.length*Vector3(0.,0.,1.)
+	V=2*(4./3)*math.pi*radius**3
+	geomInert=(2./5.)*V*radius**2+b.shape.length*b.shape.length*2*(4./3)*math.pi*radius**3
+	b.state=ChainedState()
+	b.state.addToChain(b.id)
+	#b.state.blockedDOFs=['rx','ry','rz']
+	_commonBodySetup(b,V,Vector3(geomInert,geomInert,geomInert),material,resetState=False)
+	b.state.pos=b.state.refPos=begin
+	b.dynamic=dynamic
+	b.mask=mask
+	b.bound=Aabb(diffuseColor=[0,1,0])
+	b.state.ori=b.state.ori.setFromTwoVectors(Vector3(0.,0.,1.),segment)
+	return b
+
 def wall(position,axis,sense=0,color=None,material=-1,mask=1):
 	"""Return ready-made wall body.
 

=== added file 'scripts/test/chained-cylinder-spring.py'
--- scripts/test/chained-cylinder-spring.py	1970-01-01 00:00:00 +0000
+++ scripts/test/chained-cylinder-spring.py	2010-07-12 08:50:35 +0000
@@ -0,0 +1,72 @@
+#--- bruno.chareyre@xxxxxxxxxxx ---
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+# Experiment beam-like behaviour with chained cylinders + CohFrict connexions
+
+from yade import utils
+young=1.0e6
+poisson=10
+density=2.60e3 
+frictionAngle=radians(30)
+O.materials.append(CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,label='mat'))
+O.dt=9e-5
+
+O.initializers=[
+	## Create bounding boxes. They are needed to zoom the 3d view properly before we start the simulation.
+	BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Cylinder_Aabb()])
+]
+
+O.engines=[
+	ForceResetter(),
+	BoundDispatcher([
+		Bo1_Cylinder_Aabb(),
+		Bo1_Sphere_Aabb()
+	]),
+	InsertionSortCollider(),
+	InteractionDispatchers(
+		[Ig2_ChainedCylinder_ChainedCylinder_ScGeom(),Ig2_Sphere_ChainedCylinder_CylScGeom()],
+		[Ip2_2xCohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True,normalCohesion=1e13,shearCohesion=1e13)],
+		[Law2_ScGeom_CohFrictPhys_ElasticPlastic(momentRotationLaw=True,label='law')]
+	),
+	## Apply gravity
+	GravityEngine(gravity=[0,-9.81,0]),
+	## Motion equation
+	NewtonIntegrator(damping=0.15),
+	PeriodicPythonRunner(iterPeriod=500,command='history()'),
+	#PeriodicPythonRunner(iterPeriod=100,command='yade.qt.center()')
+]
+
+#Generate a spiral
+Ne=200
+for i in range(0, Ne):
+	omega=60.0/float(Ne); hy=0.10; hz=0.15;
+	px=float(i)*(omega/60.0); py=sin(float(i)*omega)*hy; pz=cos(float(i)*omega)*hz;
+	px2=float(i+1.)*(omega/60.0); py2=sin(float(i+1.)*omega)*hy; pz2=cos(float(i+1.)*omega)*hz;
+	O.bodies.append(utils.chCylinder(begin=Vector3(pz,py,px), radius=0.005,end=Vector3(pz2,py2,px2),color=Vector3(0.6,0.5,0.5)))
+
+def outp(id=1):
+	for i in O.interactions:
+		if i.id1 == 1:
+			print i.phys.shearForce
+			print i.phys.normalForce
+			return  i
+
+O.bodies[Ne-1].state.blockedDOFs=['x','y','z','rx','ry','rz']
+yade.qt.View();
+
+
+ #plot some results
+from math import *
+from yade import plot
+plot.plots={'t':('pos1','|||','vel1')}
+def history():
+  	plot.addData(pos1=O.bodies[0].state.pos[1], # potential elastic energy
+		     vel1=O.bodies[0].state.vel[1],
+		     t=O.time)
+
+#yade.qt.Renderer().bound=True
+#O.run(1,True)
+plot.plot()
+#O.bodies[0].state.angVel=Vector3(0.05,0,0)
+