yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2490
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Fri 2010-10-15 18:49:26 +0200
message:
  1. Add Vector6r wrapper (both c++ and python)
  2. Fix prefix bug in initialization braking debian packages
  3. Fix renames in py/test which made regression tests fail
  4. Update LawTester and L3Geom for 6 dofs (not yet tested, please do not use now)
modified:
  core/main/main.py.in
  debian/rules
  lib/base/Math.hpp
  pkg/dem/DomainLimiter.cpp
  pkg/dem/DomainLimiter.hpp
  pkg/dem/L3Geom.cpp
  pkg/dem/L3Geom.hpp
  py/__init__.py.in
  py/config.py.in
  py/mathWrap/miniEigen.cpp
  py/tests/wrapper.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 'core/main/main.py.in'
--- core/main/main.py.in	2010-10-13 14:03:59 +0000
+++ core/main/main.py.in	2010-10-15 16:49:26 +0000
@@ -211,7 +211,9 @@
 		# contrary to display.Display, _BaseDisplay does not check for extensions and that avoids spurious message "Xlib.protocol.request.QueryExtension" (bug?)
 		Xlib.display._BaseDisplay();
 		yade.runtime.hasDisplay=True
-	except Xlib.error.DisplayError:
+	except: 
+		# usually Xlib.error.DisplayError, but there can be Xlib.error.XauthError etc as well
+		# let's just pretend any exception means the display would not work
 		gui=None
 
 # run remote access things, before actually starting the user session

=== modified file 'debian/rules'
--- debian/rules	2010-10-11 16:28:49 +0000
+++ debian/rules	2010-10-15 16:49:26 +0000
@@ -48,12 +48,13 @@
 	dh_installdirs
 	## files VERSION and _VERSION contain snapshot version without/with leading '-' respectively (like svn1113 and -svn1113)
 	#debug build
+	# remove the profile, if existing, to have clean config
+	rm -f scons.profile-deb
 	# the last argument specified that we only want files in lib installed -- common files (the executable) are already in the optimized package
-	NO_SCONS_GET_RECENT= scons profile=default buildPrefix=debian runtimePREFIX=/usr version=${VERSION} brief=0 chunkSize=1 jobs=1 features=vtk,gts,opengl,openmp,qt4 PREFIX=debian/yade${_VERSION}-dbg/usr variant='' optimize=0 march= debug=1 CPPPATH=/usr/include/vtk-5.0:/usr/include/vtk-5.2:/usr/include/vtk-5.4:/usr/include/eigen2 QT4DIR=/usr/share/qt4 debian/yade${_VERSION}-dbg/usr/lib
+	NO_SCONS_GET_RECENT= scons profile=deb buildPrefix=debian runtimePREFIX=/usr version=${VERSION} brief=0 chunkSize=1 jobs=1 features=vtk,gts,opengl,openmp,qt4 PREFIX=debian/yade${_VERSION}-dbg/usr variant='' optimize=0 march= debug=1 CPPPATH=/usr/include/vtk-5.0:/usr/include/vtk-5.2:/usr/include/vtk-5.4:/usr/include/eigen2 QT4DIR=/usr/share/qt4 debian/yade${_VERSION}-dbg/usr/lib
 	#optimized build
-	NO_SCONS_GET_RECENT= scons profile=default PREFIX=debian/yade${_VERSION}/usr variant='' optimize=1 debug=0
+	NO_SCONS_GET_RECENT= scons profile=deb PREFIX=debian/yade${_VERSION}/usr variant='' optimize=1 debug=0
 	#install platform-independent files (docs, scripts, examples)
-	#NO_SCONS_GET_RECENT= scons profile=deb PREFIX=debian/yade${_VERSION}/usr debian/yade${_VERSION}/usr/share/doc/yade${_VERSION}-doc
 	mkdir -p debian/yade${_VERSION}/usr/share/doc/yade${_VERSION}
 	cp -r examples scripts debian/yade${_VERSION}/usr/share/doc/yade${_VERSION}
 	# copy pixmap files in relevant directories
@@ -70,7 +71,8 @@
 	dh_testdir
 	dh_testroot
 	YADE_PREFIX=debian/yade${_VERSION}/usr debian/yade${_VERSION}/usr/bin/yade${_VERSION} --test
-	YADE_PREFIX=debian/yade${_VERSION}/usr debian/yade${_VERSION}/usr/bin/yade${_VERSION} --debug --test
+	# prefix is for libs, but the executable is only in the optimized package directory
+	YADE_PREFIX=debian/yade${_VERSION}-dbg/usr debian/yade${_VERSION}/usr/bin/yade${_VERSION} --debug --test
 
 # Build architecture-independent files here.
 binary-indep: build install

=== modified file 'lib/base/Math.hpp'
--- lib/base/Math.hpp	2010-09-10 11:41:37 +0000
+++ lib/base/Math.hpp	2010-10-15 16:49:26 +0000
@@ -40,6 +40,7 @@
 // templates of those types with single parameter are not possible (for compat with Wm3), use macros for now
 #define VECTOR2_TEMPLATE(Scalar) Eigen::Matrix<Scalar,2,1>
 #define VECTOR3_TEMPLATE(Scalar) Eigen::Matrix<Scalar,3,1>
+#define VECTOR6_TEMPLATE(Scalar) Eigen::Matrix<Scalar,6,1>
 #define MATRIX3_TEMPLATE(Scalar) Eigen::Matrix<Scalar,3,3>
 // this would be the proper way, but only works in c++-0x (not yet supported by gcc (4.5))
 #if 0
@@ -55,6 +56,7 @@
 typedef VECTOR2_TEMPLATE(Real) Vector2r;
 typedef VECTOR3_TEMPLATE(int) Vector3i;
 typedef VECTOR3_TEMPLATE(Real) Vector3r;
+typedef VECTOR6_TEMPLATE(Real) Vector6r;
 typedef MATRIX3_TEMPLATE(Real) Matrix3r;
 
 typedef Eigen::Quaternion<Real> Quaternionr;
@@ -64,6 +66,7 @@
 // io
 template<class Scalar> std::ostream & operator<<(std::ostream &os, const VECTOR2_TEMPLATE(Scalar)& v){ os << v.x()<<" "<<v.y(); return os; };
 template<class Scalar> std::ostream & operator<<(std::ostream &os, const VECTOR3_TEMPLATE(Scalar)& v){ os << v.x()<<" "<<v.y()<<" "<<v.z(); return os; };
+template<class Scalar> std::ostream & operator<<(std::ostream &os, const VECTOR6_TEMPLATE(Scalar)& v){ os << v[0]<<" "<<v[1]<<" "<<v[2]<<" "<<v[3]<<" "<<v[4]<<" "<<v[5]; return os; };
 template<class Scalar> std::ostream & operator<<(std::ostream &os, const Eigen::Quaternion<Scalar>& q){ os<<q.w()<<" "<<q.x()<<" "<<q.y()<<" "<<q.z(); return os; };
 // operators
 //template<class Scalar> VECTOR3_TEMPLATE(Scalar) operator*(Scalar s, const VECTOR3_TEMPLATE(Scalar)& v) {return v*s;}
@@ -76,6 +79,8 @@
 template<typename Scalar> bool operator!=(const Quaternion<Scalar>& u, const Quaternion<Scalar>& v){ return !(u==v); }
 template<typename Scalar> bool operator==(const MATRIX3_TEMPLATE(Scalar)& m, const MATRIX3_TEMPLATE(Scalar)& n){ for(int i=0;i<3;i++)for(int j=0;j<3;j++)if(m(i,j)!=n(i,j)) return false; return true; }
 template<typename Scalar> bool operator!=(const MATRIX3_TEMPLATE(Scalar)& m, const MATRIX3_TEMPLATE(Scalar)& n){ return !(m==n); }
+template<typename Scalar> bool operator==(const VECTOR6_TEMPLATE(Scalar)& u, const VECTOR6_TEMPLATE(Scalar)& v){ return u[0]==v[0] && u[1]==v[1] && u[2]==v[2] && u[3]==v[3] && u[4]==v[4] && u[5]==v[5]; }
+template<typename Scalar> bool operator!=(const VECTOR6_TEMPLATE(Scalar)& u, const VECTOR6_TEMPLATE(Scalar)& v){ return !(u==v); }
 template<typename Scalar> bool operator==(const VECTOR3_TEMPLATE(Scalar)& u, const VECTOR3_TEMPLATE(Scalar)& v){ return u.x()==v.x() && u.y()==v.y() && u.z()==v.z(); }
 template<typename Scalar> bool operator!=(const VECTOR3_TEMPLATE(Scalar)& u, const VECTOR3_TEMPLATE(Scalar)& v){ return !(u==v); }
 template<typename Scalar> bool operator==(const VECTOR2_TEMPLATE(Scalar)& u, const VECTOR2_TEMPLATE(Scalar)& v){ return u.x()==v.x() && u.y()==v.y(); }
@@ -180,6 +185,7 @@
 BOOST_IS_BITWISE_SERIALIZABLE(Vector2i);
 BOOST_IS_BITWISE_SERIALIZABLE(Vector3r);
 BOOST_IS_BITWISE_SERIALIZABLE(Vector3i);
+BOOST_IS_BITWISE_SERIALIZABLE(Vector6r);
 BOOST_IS_BITWISE_SERIALIZABLE(Quaternionr);
 BOOST_IS_BITWISE_SERIALIZABLE(Se3r);
 BOOST_IS_BITWISE_SERIALIZABLE(Matrix3r);
@@ -213,6 +219,12 @@
 }
 
 template<class Archive>
+void serialize(Archive & ar, Vector6r & g, const unsigned int version){
+	Real &v0=g[0], &v1=g[1], &v2=g[2], &v3=g[3], &v4=g[4], &v5=g[5];
+	ar & BOOST_SERIALIZATION_NVP(v0) & BOOST_SERIALIZATION_NVP(v1) & BOOST_SERIALIZATION_NVP(v2) & BOOST_SERIALIZATION_NVP(v3) & BOOST_SERIALIZATION_NVP(v4) & BOOST_SERIALIZATION_NVP(v5);
+}
+
+template<class Archive>
 void serialize(Archive & ar, Quaternionr & g, const unsigned int version)
 {
 	Real &w=g.w(), &x=g.x(), &y=g.y(), &z=g.z();

=== modified file 'pkg/dem/DomainLimiter.cpp'
--- pkg/dem/DomainLimiter.cpp	2010-10-14 17:21:34 +0000
+++ pkg/dem/DomainLimiter.cpp	2010-10-15 16:49:26 +0000
@@ -31,16 +31,18 @@
 
 void LawTester::postLoad(LawTester&){
 	if(ids.size()!=2) throw std::invalid_argument("LawTester.ids: exactly two values must be given.");
-	if(path.empty()) throw invalid_argument("LawTester.path: at least one point must be given.");
+	if(path.empty() && rotPath.empty()) throw invalid_argument("LawTester.{path,rotPath}: at least one point must be given.");
 	if(pathSteps.empty()) throw invalid_argument("LawTester.pathSteps: at least one value must be given.");
+	size_t pathSize=max(path.size(),rotPath.size());
 	// update path points
-	_pathV.clear(); _pathV.push_back(Vector3r::Zero());
-	for(size_t i=0; i<path.size(); i++) _pathV.push_back(path[i]);
+	_pathU.clear(); _pathU.push_back(Vector3r::Zero());
+	_pathR.clear(); _pathR.push_back(Vector3r::Zero());
+	for(size_t i=0; i<pathSize; i++) { _pathU.push_back(i<path.size()?path[i]:*(path.rbegin())); _pathR.push_back(i<rotPath.size()?rotPath[i]:*(rotPath.rbegin())); }
 	// update time points from distances, repeat last distance if shorter than path
 	_pathT.clear(); _pathT.push_back(0);
 	for(size_t i=0; i<pathSteps.size(); i++) _pathT.push_back(_pathT[i]+pathSteps[i]);
 	int lastDist=pathSteps[pathSteps.size()-1];
-	for(size_t i=pathSteps.size(); i<path.size(); i++) _pathT.push_back(*(_pathT.rbegin())+lastDist);
+	for(size_t i=pathSteps.size(); i<pathSize; i++) _pathT.push_back(*(_pathT.rbegin())+lastDist);
 }
 
 void LawTester::action(){
@@ -63,9 +65,10 @@
 	ScGeom* scGeom=dynamic_cast<ScGeom*>(I->geom.get());
 	Dem3DofGeom* d3dGeom=dynamic_cast<Dem3DofGeom*>(I->geom.get());
 	L3Geom* l3Geom=dynamic_cast<L3Geom*>(I->geom.get());
+	L6Geom* l6Geom=dynamic_cast<L6Geom*>(I->geom.get());
 	//NormShearPhys* phys=dynamic_cast<NormShearPhys*>(I->phys.get());			//Disabled because of warning
 	if(!gsc) throw std::invalid_argument("LawTester: IGeom of "+strIds+" not a GenericSpheresContact.");
-	if(!scGeom && !d3dGeom && !l3Geom) throw std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom, nor L3Geom.");
+	if(!scGeom && !d3dGeom && !l3Geom) throw std::invalid_argument("LawTester: IGeom of "+strIds+" is neither ScGeom, nor Dem3DofGeom, nor L3Geom (or L6Geom).");
 	assert(!((bool)scGeom && (bool)d3dGeom && (bool)l3Geom)); // nonsense
 	// get body objects
 	State *state1=Body::byId(id1,scene)->state.get(), *state2=Body::byId(id2,scene)->state.get();
@@ -97,6 +100,7 @@
 				scGeom->rotate(shearTot);
 				shearTot+=scGeom->shearIncrement();
 				ptGeom=Vector3r(-scGeom->penetrationDepth,shearTot.dot(axY),shearTot.dot(axZ));
+				rotGeom=Vector3r(NaN,NaN,NaN);
 			}
 			else{ // d3dGeom
 				throw runtime_error("LawTester: Dem3DofGeom not yet supported.");
@@ -111,6 +115,7 @@
 		trsf=l3Geom->trsf;
 		axX=trsf.col(0); axY=trsf.col(1); axZ=trsf.col(2);
 		ptGeom=l3Geom->u;
+		if(l6Geom) rotGeom=l6Geom->phi;
 	}
 	trsfQ=Quaternionr(trsf);
 	contPt=gsc->contactPoint;
@@ -118,14 +123,17 @@
 	renderLength=.5*refLength;
 	
 	// here we go ahead, finally
-	Vector3r val=linearInterpolate<Vector3r,int>(step,_pathT,_pathV,_interpPos);
-	Vector3r valDiff=val-prevVal; prevVal=val;
+	Vector3r u=linearInterpolate<Vector3r,int>(step,_pathT,_pathU,_interpPos);
+	Vector3r phi=linearInterpolate<Vector3r,int>(step,_pathT,_pathR,_interpPosRot);
+	Vector3r dU=u-uPrev; uPrev=u;
+	Vector3r dPhi=phi-phiPrev; phiPrev=phi;
 	if(displIsRel){
-		LOG_DEBUG("Relative diff is "<<valDiff<<" (will be normalized by "<<gsc->refR1+gsc->refR2<<")");
-		valDiff*=refLength;
+		LOG_DEBUG("Relative displacement diff is "<<dU<<" (will be normalized by "<<gsc->refR1+gsc->refR2<<")");
+		dU*=refLength;
 	}
-	LOG_DEBUG("Absolute diff is "<<valDiff);
-	ptOurs+=valDiff;
+	LOG_DEBUG("Absolute diff is: displacement "<<dU<<", rotation "<<dPhi);
+	ptOurs+=dU;
+	rotOurs+=dPhi;
 
 	// reset velocities where displacement is controlled
 	//for(int i=0; i<3; i++){ if(forceControl[i]==0){ state1.vel[i]=0; state2.vel[i]=0; }
@@ -137,29 +145,31 @@
 		int sign=(i==0?-1:1);
 		Real weight=(i==0?1-idWeight:idWeight);
 		Real radius=(i==0?gsc->refR1:gsc->refR2);
-		// signed and weighted displacement to be applied on this sphere (reversed for #0)
+		// signed and weighted displacement/rotation to be applied on this sphere (reversed for #0)
 		// some rotations must cancel the sign, by multiplying by sign again
-		Vector3r diff=sign*valDiff*weight;
+		Vector3r ddU=sign*dU*weight;
+		Vector3r ddPhi=sign*dPhi*weight;
+		vel[i]=angVel[i]=Vector3r::Zero();
 
 		// normal displacement
 
-		vel[i]=axX*diff[0]/scene->dt;
+		vel[i]+=axX*ddU[0]/scene->dt;
 
 		// shear rotation
 
-		//   multiplication by sign cancels sign in diff, since rotation is non-symmetric (to increase shear, both spheres have the same rotation)
+		//   multiplication by sign cancels sign in ddU, since rotation is non-symmetric (to increase shear, both spheres have the same rotation)
 		//   (unlike shear displacement, which is symmetric)
 		// rotation around Z (which gives y-shear) must be inverted: +ry gives +εzm while -rz gives +εy
-		Real rotZ=-sign*rotWeight*diff[1]/radius, rotY=sign*rotWeight*diff[2]/radius;
-		angVel[i]=(rotY*axY+rotZ*axZ)/scene->dt;
+		Real rotZ=-sign*rotWeight*ddU[1]/radius, rotY=sign*rotWeight*ddU[2]/radius;
+		angVel[i]+=(rotY*axY+rotZ*axZ)/scene->dt;
 
 		// shear displacement
 
-		// angle that is traversed by a sphere in order to give desired diff when displaced on the branch of r1+r2
+		// angle that is traversed by a sphere in order to give desired ddU when displaced on the branch of r1+r2
 		// FIXME: is the branch value correct here?!
-		Real arcAngleY=atan((1-rotWeight)*diff[1]/radius), arcAngleZ=atan((1-rotWeight)*diff[2]/radius);
+		Real arcAngleY=atan((1-rotWeight)*ddU[1]/radius), arcAngleZ=atan((1-rotWeight)*ddU[2]/radius);
 		// same, but without the atan, which can be disregarded for small increments:
-		//    Real arcAngleY=(1-rotWeight)*diff[1]/radius, arcAngleZ=(1-rotWeight)*diff[2]/radius; 
+		//    Real arcAngleY=(1-rotWeight)*ddU[1]/radius, arcAngleZ=(1-rotWeight)*ddU[2]/radius; 
 		vel[i]+=axY*radius*sin(arcAngleY)/scene->dt; vel[i]+=axZ*radius*sin(arcAngleZ)/scene->dt;
 
 		// compensate distance increase caused by motion along the perpendicular axis
@@ -167,10 +177,11 @@
 		// and the compensation is always in the -εx sense (-sign → +1 for #0, -1 for #1)
 		vel[i]+=-sign*axX*radius*((1-cos(arcAngleY))+(1-cos(arcAngleZ)))/scene->dt;
 
+		// rotation
+		angVel[i]+=trsf*ddPhi;
+
 		LOG_DEBUG("vel="<<vel[i]<<", angVel="<<angVel[i]<<", rotY,rotZ="<<rotY<<","<<rotZ<<", arcAngle="<<arcAngleY<<","<<arcAngleZ<<", sign="<<sign<<", weight="<<weight);
 
-		//Vector3r vel=(rotY*axY+rotZ*axZ)/scene->dt;
-		//Vector3r linVel=axX*valDiff[0]/scene->dt;
 	}
 	state1->vel=vel[0]; state1->angVel=angVel[0];
 	state2->vel=vel[1]; state2->angVel=angVel[1];
@@ -228,7 +239,7 @@
 	glTranslatev(Vector3r(0,tester->ptOurs[1],tester->ptOurs[2]));
 
 
-	const int t(tester->step); const vector<int>& TT(tester->_pathT); const vector<Vector3r>& VV(tester->_pathV);
+	const int t(tester->step); const vector<int>& TT(tester->_pathT); const vector<Vector3r>& VV(tester->_pathU);
 	size_t numSegments=TT.size();
 	const Vector3r colorBefore=Vector3r(.7,1,.7), colorAfter=Vector3r(1,.7,.7);
 

=== modified file 'pkg/dem/DomainLimiter.hpp'
--- pkg/dem/DomainLimiter.hpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/DomainLimiter.hpp	2010-10-15 16:49:26 +0000
@@ -21,23 +21,30 @@
 		void postLoad(LawTester&);
 	DECLARE_LOGGER;
 	YADE_CLASS_BASE_DOC_ATTRS(LawTester,PartialEngine,"Prescribe and apply deformations of an interaction in terms of normal and shear displacements. See :ysrc:`scripts/test/law-test.py`. ",
-		((vector<Vector3r>,path,,Attr::triggerPostLoad,"Loading path, where each Vector3 contains desired normal displacement and two components of the shear displacement (in local coordinate system, which is being tracked automatically."))
+		((vector<Vector3r>,path,,Attr::triggerPostLoad,"Loading path, where each Vector3 contains desired normal displacement and two components of the shear displacement (in local coordinate system, which is being tracked automatically. If shorter than :yref:`rotPath<LawTester.rotPath>`, the last value is repeated."))
+		((vector<Vector3r>,rotPath,,Attr::triggerPostLoad,"Rotational components of the loading path, where each item contains torsion and two bending rotations in local coordinates. If shorter than :yref:`path<LawTester.path>`, the last value is repeated."))
 		((vector<string>,hooks,,,"Python commands to be run when the corresponding point in path is reached, before doing other things in that particular step. See also :yref:`doneHook<LawTester.doneHook>`. "))
-		((Vector3r,ptOurs,Vector3r::Zero(),,"Current absolute configuration -- computed by ourselves from applied increments; should correspond to posTot."))
-		((Vector3r,shearTot,Vector3r::Zero(),,"Current total shear, in global coordinates (computed differently for :yref:`ScGeom` and :yref:`Dem3DofGeom`)."))
-		((Vector3r,ptGeom,Vector3r::Zero(),,"Current absolute configuration, in local coordinates; computed from *shearTot*, indirectly  from :yref:`IGeom`. |yupdate|"))
+		((Vector3r,ptOurs,Vector3r::Zero(),,"Current displacement, computed by ourselves from applied increments; should correspond to ptGeom."))
+		((Vector3r,ptGeom,Vector3r::Zero(),,"Current displacement, as computed by the geometry functor"))
+		((Vector3r,shearTot,Vector3r::Zero(),,"Current displacement in global coordinates."))
+		((Vector3r,rotOurs,Vector3r::Zero(),,"Current rotation, computed by ourselves from applied increments; should correspond to rotGeom."))
+		((Vector3r,rotGeom,Vector3r::Zero(),,"Current rotation, as computed by the geometry functor"))
+		((Vector3r,rotTot,Vector3r::Zero(),,"Current rotation in global coordinates."))
 		((bool,displIsRel,true,,"Whether displacement values in *path* are normalized by reference contact length (r1+r2 for 2 spheres)."))
 		((Vector3i,forceControl,Vector3i::Zero(),,"Select which components of path (non-zero value) have force (stress) rather than displacement (strain) meaning."))
 		((vector<int>,pathSteps,((void)"(constant step)",vector<int>(1,1)),Attr::triggerPostLoad,"Step number for corresponding values in :yref:`path<LawTester.path>`; if shorter than path, distance between last 2 values is used for the rest."))
 		((vector<int>,_pathT,,(Attr::readonly|Attr::noSave),"Time value corresponding to points on path, computed from *pathSteps*. Length is always the same as path."))
-		((vector<Vector3r>,_pathV,,(Attr::readonly|Attr::noSave),"Path values, computed from *path* by appending zero initial value."))
+		((vector<Vector3r>,_pathU,,(Attr::readonly|Attr::noSave),"Displacement path values, computed from *path* by appending zero initial value (and possibly repeating the last value to match the length of rotation path."))
+		((vector<Vector3r>,_pathR,,(Attr::readonly|Attr::noSave),"Rotation path values, computed from *path* by appending zero initial value (and possibly repeating the last value to match the length of displacement path)."))
 		((shared_ptr<Interaction>,I,,(Attr::hidden),"Interaction object being tracked."))
 		((Vector3r,axX,,,"Local x-axis in global coordinates (normal of the contact) |yupdate|"))
 		((Vector3r,axY,,,"Local y-axis in global coordinates; perpendicular to axX; initialized arbitrarily, but tracked to be consistent. |yupdate|"))
 		((Vector3r,axZ,,Attr::noSave,"Local z-axis in global coordinates; computed from axX and axY. |yupdate|"))
 		((Matrix3r,trsf,,Attr::noSave,"Transformation matrix for the local coordinate system. |yupdate|"))
-		((size_t,_interpPos,0,(Attr::readonly),"Position for the interpolation routine."))
-		((Vector3r,prevVal,Vector3r::Zero(),(Attr::readonly),"Value reached in the previous step."))
+		((size_t,_interpPos,0,(Attr::readonly|Attr::hidden),"State cookie for the interpolation routine (displacements)."))
+		((size_t,_interpPosRot,0,(Attr::readonly|Attr::hidden),"State cookie for the interpolation routine (rotations)."))
+		((Vector3r,uPrev,Vector3r::Zero(),(Attr::readonly),"Displacement value reached in the previous step."))
+		((Vector3r,phiPrev,Vector3r::Zero(),(Attr::readonly),"Rotation value reached in the previous step."))
 		((Quaternionr,trsfQ,,Attr::noSave,"Transformation quaterion for the local coordinate system. |yupdate|"))
 		((int,step,0,,"Step number in which this engine is active; determines position in path, using pathSteps."))
 		((string,doneHook,,,"Python command (as string) to run when end of the path is achieved. If empty, the engine will be set :yref:`dead<Engine.dead>`."))

=== modified file 'pkg/dem/L3Geom.cpp'
--- pkg/dem/L3Geom.cpp	2010-10-14 17:21:34 +0000
+++ pkg/dem/L3Geom.cpp	2010-10-15 16:49:26 +0000
@@ -2,11 +2,22 @@
 #include<yade/pkg-dem/L3Geom.hpp>
 #include<yade/pkg-common/Sphere.hpp>
 
-YADE_PLUGIN((L3Geom)(Ig2_Sphere_Sphere_L3Geom_Inc)(Law2_L3Geom_FrictPhys_Linear));
+YADE_PLUGIN((L3Geom)(Ig2_Sphere_Sphere_L3Geom_Inc)(Law2_L3Geom_FrictPhys_Linear)(Law2_L6Geom_FrictPhys_Linear));
 
 L3Geom::~L3Geom(){}
+L6Geom::~L6Geom(){}
 
 bool Ig2_Sphere_Sphere_L3Geom_Inc::go(const shared_ptr<Shape>& s1, const shared_ptr<Shape>& s2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& I){
+	return genericGo(/*is6Dof*/false,s1,s2,state1,state2,shift2,force,I);
+};
+
+bool Ig2_Sphere_Sphere_L6Geom_Inc::go(const shared_ptr<Shape>& s1, const shared_ptr<Shape>& s2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& I){
+	return genericGo(/*is6Dof*/true,s1,s2,state1,state2,shift2,force,I);
+};
+
+
+bool Ig2_Sphere_Sphere_L3Geom_Inc::genericGo(bool is6Dof, const shared_ptr<Shape>& s1, const shared_ptr<Shape>& s2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& I){
+
 	const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
 	const Real& r1(static_pointer_cast<Sphere>(s1)->radius); const Real& r2(static_pointer_cast<Sphere>(s2)->radius);
 	Vector3r relPos=(se32.position+shift2)-se31.position;
@@ -21,8 +32,10 @@
 	Vector3r contPt=se31.position+(r1+0.5*uN)*normal;
 
 	// create geometry
-	if(!I->geom){ 
-		I->geom=shared_ptr<L3Geom>(new L3Geom); const shared_ptr<L3Geom>& g(static_pointer_cast<L3Geom>(I->geom));
+	if(!I->geom){
+		if(is6Dof) I->geom=shared_ptr<L6Geom>(new L6Geom);
+		else       I->geom=shared_ptr<L3Geom>(new L3Geom);
+		const shared_ptr<L3Geom>& g(static_pointer_cast<L3Geom>(I->geom));
 		g->contactPoint=contPt;
 		g->refR1=r1; g->refR2=r1;
 		g->normal=normal; const Vector3r& locX(g->normal);
@@ -31,6 +44,7 @@
 		Vector3r locZ=g->normal.cross(locY);
 		g->trsf.col(0)=locX; g->trsf.col(1)=locY; g->trsf.col(2)=locZ;
 		g->u=Vector3r(uN,0,0); // zero shear displacement
+		// L6Geom::phi is initialized to Vector3r::Zero() automatically
 		return true;
 	}
 	
@@ -47,7 +61,8 @@
 	// contrary to what ScGeom::precompute does now (r2486), we take average normal, i.e. .5*(prevNormal+currNormal),
 	// so that all terms in the equation are in the previous mid-step
 	// the re-normalization might not be necessary for very small increments, but better do it
-	Vector3r avgNormal=.5*(prevNormal+currNormal); avgNormal.normalize();
+	Vector3r avgNormal=(approxMask|APPROX_NO_MID_NORMAL) ? prevNormal : .5*(prevNormal+currNormal);
+	if(!(approxMask|APPROX_NO_RENORM_MID_NORMAL) && !(approxMask|APPROX_NO_MID_NORMAL)) avgNormal.normalize(); // normalize only if used and if requested via approxMask
 	// twist vector of the normal from the last step
 	Vector3r normTwistVec=avgNormal*scene->dt*.5*avgNormal.dot(state1.angVel+state2.angVel);
 	// compute relative velocity
@@ -76,6 +91,7 @@
 	for(int i=1; i<3; i++){
 		currTrsf.col(i)=prevTrsf.col(i)-prevTrsf.col(i).cross(normRotVec)-prevTrsf.col(i).cross(normTwistVec);
 	}
+	if(!(approxMask | APPROX_NO_RENORM_TRSF)){ /* renormalizing quternion is faster*/ currTrsf=Matrix3r(Quaternionr(currTrsf).normalized()); }
 
 	/* Previous local trsf u'⁻ must be updated to current u'⁰. We have transformation T⁻ and T⁰,
 		δ(a) denotes increment of a as defined above.  Two possibilities:
@@ -88,7 +104,8 @@
 			This could be perhaps simplified by using T⁰ or T⁻ since they will not differ much,
 			but it would have to be verified somehow.
 	*/
-	Quaternionr midTrsf=Quaternionr(prevTrsf).slerp(.5,Quaternionr(currTrsf));
+	// if requested via approxMask, just use prevTrsf
+	Quaternionr midTrsf=(approxMask|APPROX_NO_MID_TRSF) ? Quaternionr(prevTrsf) : Quaternionr(prevTrsf).slerp(.5,Quaternionr(currTrsf));
 	
 	// updates of geom here
 
@@ -102,6 +119,13 @@
 	g->normal=currNormal;
 	g->contactPoint=contPt;
 
+	if(is6Dof){
+		const shared_ptr<L6Geom> g6=static_pointer_cast<L6Geom>(g);
+		// update phi, from the difference of angular velocities
+		// the difference is transformed to local coord using the midTrsf transformation
+		g6->phi+=midTrsf*(scene->dt*(state2.angVel-state1.angVel));
+	}
+
 	return true;
 };
 
@@ -118,3 +142,14 @@
 	applyForceAtBranch(globalF,I->getId1(),geom->normal*(geom->refR1+.5*geom->u[0]),I->getId2(),-geom->normal*(geom->refR1+.5*geom->u[0]));
 }
 
+void Law2_L6Geom_FrictPhys_Linear::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
+	// same as for L3Geom, except of the torque term
+	const shared_ptr<L6Geom> geom=static_pointer_cast<L6Geom>(ig); const shared_ptr<FrictPhys> phys=static_pointer_cast<FrictPhys>(ip);
+	Vector3r localF=geom->u.cwise()*Vector3r(phys->kn,phys->ks,phys->ks);
+	Vector3r localT=charLen*(geom->phi.cwise()*Vector3r(phys->kn,phys->ks,phys->ks));
+	Quaternionr invTrsf=Quaternionr(geom->trsf).conjugate();
+	Vector3r globalF=invTrsf*localF, globalT=invTrsf*localT;
+	phys->normalForce=geom->normal*globalF.dot(geom->normal); phys->shearForce=globalF-phys->normalForce;
+	applyForceAtBranch(globalF,I->getId1(),geom->normal*(geom->refR1+.5*geom->u[0]),I->getId2(),-geom->normal*(geom->refR1+.5*geom->u[0]));
+	scene->forces.addTorque(I->getId1(),globalT); scene->forces.addTorque(I->getId2(),-globalT);
+}

=== modified file 'pkg/dem/L3Geom.hpp'
--- pkg/dem/L3Geom.hpp	2010-10-14 17:21:34 +0000
+++ pkg/dem/L3Geom.hpp	2010-10-15 16:49:26 +0000
@@ -8,8 +8,8 @@
 #include<yade/pkg-dem/FrictPhys.hpp>
 
 struct L3Geom: public GenericSpheresContact{
-	Real& uN;
-	Vector2r& uT; 
+	const Real& uN;
+	const Vector2r& uT; 
 	virtual ~L3Geom();
 	YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(L3Geom,GenericSpheresContact,"Geometry of contact given in local coordinates with 3 degress of freedom: normal and two in shear plane. [experimental]",
 		((Vector3r,u,Vector3r::Zero(),,"Displacement components, in local coordinates. |yupdate|"))
@@ -22,7 +22,7 @@
 		((Matrix3r,trsf,Matrix3r::Identity(),,"Transformation (rotation) from global to local coordinates. (the translation part is in :yref:`GenericSpheresContact.contactPoint`)"))
 		,
 		/*init*/
-		((uN,u[0])) ((uT,*((Vector2r*)&u[1])))
+		((uN,u[0])) ((uT,Vector2r::Map(&u[1])))
 		,
 		/*ctor*/ createIndex();
 		, /*py*/
@@ -33,29 +33,55 @@
 };
 REGISTER_SERIALIZABLE(L3Geom);
 
-/*
-class L3Geom_Sphere_Sphere: public L3Geom{
-	YADE_CLASS_BASE_DOC_ATTRS_CTOR(L3Geom_Sphere_Sphere,L3Geom,"Local geometry of contact between two spherical particles with 3 DoFs. [experimental]",
-		((
-
-}
-*/
+
+struct L6Geom: public L3Geom{
+	virtual ~L6Geom();
+	YADE_CLASS_BASE_DOC_ATTRS(L6Geom,L3Geom,"Geoemtric of contact in local coordinates with 6 degrees of freedom. [experimental]",
+		((Vector3r,phi,Vector3r::Zero(),,"Rotation components, in local coordinates. |yupdate|"))
+	);
+};
+REGISTER_SERIALIZABLE(L6Geom);
+
 
 struct Ig2_Sphere_Sphere_L3Geom_Inc: public IGeomFunctor{
 		virtual bool go(const shared_ptr<Shape>& s1, const shared_ptr<Shape>& s2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& I);
+		virtual bool genericGo(bool is6Dof, const shared_ptr<Shape>& s1, const shared_ptr<Shape>& s2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& I);
+
+	enum { APPROX_NO_RENORM_TRSF=1, APPROX_NO_MID_TRSF=2, APPROX_NO_MID_NORMAL=4, APPROX_NO_RENORM_MID_NORMAL=8 };
 
 	YADE_CLASS_BASE_DOC_ATTRS(Ig2_Sphere_Sphere_L3Geom_Inc,IGeomFunctor,"Incrementally compute :yref:`L3Geom` for contact of 2 spheres.\n\n.. note:: The initial value of *u[0]* (normal displacement) might be non-zero, with or without *distFactor*, since it is given purely by sphere's geometry. If you want to set \"equilibrium distance\", do it in the contact law as explained in :yref:`L3Geom.u0`.",
 		((bool,noRatch,true,,"See :yref:`ScGeom.avoidGranularRatcheting`."))
 		((Real,distFactor,1,,"Create interaction if spheres are not futher than distFactor*(r1+r2)."))
+		((int,approxMask,0,,"Selectively enable geometrical approximations (bitmask); add the values for approximations to be enabled.\n\n1: do not renormalize transformation matrix at every step\n2: use previous transformation to transform velocities (which are known at mid-steps), instead of mid-step transformation computed as quaternion slerp at t=0.5.\n4: do not take average (mid-step) normal when computing relative shear displacement, use previous value instead\n8: do not re-normalize average (mid-step) normal, if used.…\nBy default, the mask is zero and neither of these approximations is used."))
 	);
 	FUNCTOR2D(Sphere,Sphere);
 	DEFINE_FUNCTOR_ORDER_2D(Sphere,Sphere);
 };
 REGISTER_SERIALIZABLE(Ig2_Sphere_Sphere_L3Geom_Inc);
 
+
+struct Ig2_Sphere_Sphere_L6Geom_Inc: public Ig2_Sphere_Sphere_L3Geom_Inc{
+	virtual bool go(const shared_ptr<Shape>& s1, const shared_ptr<Shape>& s2, const State& state1, const State& state2, const Vector3r& shift2, const bool& force, const shared_ptr<Interaction>& I);
+	YADE_CLASS_BASE_DOC(Ig2_Sphere_Sphere_L6Geom_Inc,Ig2_Sphere_Sphere_L3Geom_Inc,"Incrementally compute :yref:`L6Geom` for contact of 2 spheres.");
+	FUNCTOR2D(Sphere,Sphere);
+	DEFINE_FUNCTOR_ORDER_2D(Sphere,Sphere);
+};
+REGISTER_SERIALIZABLE(Ig2_Sphere_Sphere_L6Geom_Inc);
+
+
 struct Law2_L3Geom_FrictPhys_Linear: public LawFunctor{
 	virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
 	FUNCTOR2D(L3Geom,FrictPhys);
 	YADE_CLASS_BASE_DOC(Law2_L3Geom_FrictPhys_Linear,LawFunctor,"Basic law for testing :yref:`L3Geom` -- linear in both normal and shear sense, without slip or breakage.");
 };
 REGISTER_SERIALIZABLE(Law2_L3Geom_FrictPhys_Linear);
+
+
+struct Law2_L6Geom_FrictPhys_Linear: public LawFunctor{
+	virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+	FUNCTOR2D(L6Geom,FrictPhys);
+	YADE_CLASS_BASE_DOC_ATTRS(Law2_L6Geom_FrictPhys_Linear,LawFunctor,"Basic law for testing :yref:`L3Geom` -- linear in both normal and shear sense, without slip or breakage.",
+		((Real,charLen,1,,"Characteristic length with the meaning of the stiffness ratios bending/shear and torsion/normal."))
+	);
+};
+REGISTER_SERIALIZABLE(Law2_L6Geom_FrictPhys_Linear);

=== modified file 'py/__init__.py.in'
--- py/__init__.py.in	2010-10-14 08:18:42 +0000
+++ py/__init__.py.in	2010-10-15 16:49:26 +0000
@@ -42,11 +42,11 @@
 plugins=[]
 # might add followlinks=True to os.walk, for python>=2.6
 for root,dirs,files in os.walk(config.libDir):
-	# for some reason, does not really work?
-	#dirs=[d for d in dirs if d!='dbg']
+	# skip any directory named 'dbg' in the non-debug build
+	if 'dbg' in dirs: dirs.remove('dbg')
 	for f in files:
 		# ouch, ugly!
-		if not config.debug and '/dbg/' in root: continue
+		#if not config.debug and '/dbg/' in root: continue
 		if not (f.startswith('lib') and f.endswith('.so')): continue
 		plugin=os.path.join(config.libDir,root,f)
 		plugins.append(plugin)

=== modified file 'py/config.py.in'
--- py/config.py.in	2010-10-11 16:28:49 +0000
+++ py/config.py.in	2010-10-15 16:49:26 +0000
@@ -8,7 +8,7 @@
 prefix='${runtimePREFIX}' if not os.environ.has_key('YADE_PREFIX') else os.environ['YADE_PREFIX']
 suffix='${SUFFIX}'
 debug=bool(${debug})
-libDir=os.path.abspath('$runtimePREFIX/lib/yade$SUFFIX'+('/dbg' if debug else ''))
+libDir=os.path.abspath(prefix+'/lib/yade$SUFFIX'+('/dbg' if debug else ''))
 confDir=os.environ['HOME']+'/.yade$SUFFIX'
 libstdcxx='${libstdcxx}'
 features='${features}'.split(',')

=== modified file 'py/mathWrap/miniEigen.cpp'
--- py/mathWrap/miniEigen.cpp	2010-09-30 18:00:41 +0000
+++ py/mathWrap/miniEigen.cpp	2010-10-15 16:49:26 +0000
@@ -16,6 +16,7 @@
 #define IDX_CHECK(i,MAX){ if(i<0 || i>=MAX) { PyErr_SetString(PyExc_IndexError, ("Index out of range 0.." + boost::lexical_cast<std::string>(MAX-1)).c_str()); bp::throw_error_already_set(); } }
 #define IDX2_CHECKED_TUPLE_INTS(tuple,max2,arr2) {int l=bp::len(tuple); if(l!=2) { PyErr_SetString(PyExc_IndexError,"Index must be integer or a 2-tuple"); bp::throw_error_already_set(); } for(int _i=0; _i<2; _i++) { bp::extract<int> val(tuple[_i]); if(!val.check()) throw std::runtime_error("Unable to convert "+boost::lexical_cast<std::string>(_i)+"-th index to int."); int v=val(); IDX_CHECK(v,max2[_i]); arr2[_i]=v; }  }
 
+void Vector6r_set_item(Vector6r & self, int idx, Real value){ IDX_CHECK(idx,6); self[idx]=value; }
 void Vector3r_set_item(Vector3r & self, int idx, Real value){ IDX_CHECK(idx,3); self[idx]=value; }
 void Vector3i_set_item(Vector3i & self, int idx, int  value){ IDX_CHECK(idx,3); self[idx]=value; }
 void Vector2r_set_item(Vector2r & self, int idx, Real value){ IDX_CHECK(idx,2); self[idx]=value; }
@@ -25,6 +26,7 @@
 void Matrix3r_set_item(Matrix3r & self, bp::tuple _idx, Real value){ int idx[2]; int mx[2]={3,3}; IDX2_CHECKED_TUPLE_INTS(_idx,mx,idx); self(idx[0],idx[1])=value; }
 void Matrix3r_set_item_linear(Matrix3r & self, int idx, Real value){ IDX_CHECK(idx,9); self(idx/3,idx%3)=value; }
 
+Real Vector6r_get_item(const Vector6r & self, int idx){ IDX_CHECK(idx,6); return self[idx]; }
 Real Vector3r_get_item(const Vector3r & self, int idx){ IDX_CHECK(idx,3); return self[idx]; }
 int  Vector3i_get_item(const Vector3i & self, int idx){ IDX_CHECK(idx,3); return self[idx]; }
 Real Vector2r_get_item(const Vector2r & self, int idx){ IDX_CHECK(idx,2); return self[idx]; }
@@ -34,6 +36,7 @@
 Real Matrix3r_get_item(Matrix3r & self, bp::tuple _idx){ int idx[2]; int mx[2]={3,3}; IDX2_CHECKED_TUPLE_INTS(_idx,mx,idx); return self(idx[0],idx[1]); }
 Real Matrix3r_get_item_linear(Matrix3r & self, int idx){ IDX_CHECK(idx,9); return self(idx/3,idx%3); }
 
+std::string Vector6r_str(const Vector3r & self){ return std::string("Vector6(")+boost::lexical_cast<std::string>(self[0])+","+boost::lexical_cast<std::string>(self[1])+","+boost::lexical_cast<std::string>(self[2])+", "+boost::lexical_cast<std::string>(self[3])+","+boost::lexical_cast<std::string>(self[4])+","+boost::lexical_cast<std::string>(self[5])+")";}
 std::string Vector3r_str(const Vector3r & self){ return std::string("Vector3(")+boost::lexical_cast<std::string>(self[0])+","+boost::lexical_cast<std::string>(self[1])+","+boost::lexical_cast<std::string>(self[2])+")";}
 std::string Vector3i_str(const Vector3i & self){ return std::string("Vector3i(")+boost::lexical_cast<std::string>(self[0])+","+boost::lexical_cast<std::string>(self[1])+","+boost::lexical_cast<std::string>(self[2])+")";}
 std::string Vector2r_str(const Vector2r & self){ return std::string("Vector2(")+boost::lexical_cast<std::string>(self[0])+","+boost::lexical_cast<std::string>(self[1])+")";}
@@ -41,6 +44,7 @@
 std::string Quaternionr_str(const Quaternionr & self){ AngleAxisr aa(self); return std::string("Quaternion((")+boost::lexical_cast<std::string>(aa.axis()[0])+","+boost::lexical_cast<std::string>(aa.axis()[1])+","+boost::lexical_cast<std::string>(aa.axis()[2])+"),"+boost::lexical_cast<std::string>(aa.angle())+")";}
 std::string Matrix3r_str(const Matrix3r & self){ std::ostringstream oss; oss<<"Matrix3("; for(int i=0; i<3; i++) for(int j=0; j<3; j++) oss<<self(i,j)<<((i==2 && j==2)?")":",")<<((i<2 && j==2)?" ":""); return oss.str(); }
 
+int Vector6r_len(){return 6;}
 int Vector3r_len(){return 3;}
 int Vector3i_len(){return 3;}
 int Vector2r_len(){return 2;}
@@ -51,6 +55,7 @@
 // pickling support
 struct Matrix3r_pickle: bp::pickle_suite{static bp::tuple getinitargs(const Matrix3r& x){ return bp::make_tuple(x(0,0),x(0,1),x(0,2),x(1,0),x(1,1),x(1,2),x(2,0),x(2,1),x(2,2));} };
 struct Quaternionr_pickle: bp::pickle_suite{static bp::tuple getinitargs(const Quaternionr& x){ return bp::make_tuple(x.w(),x.x(),x.y(),x.z());} };
+struct Vector6r_pickle: bp::pickle_suite{static bp::tuple getinitargs(const Vector3r& x){ return bp::make_tuple(x[0],x[1],x[2],x[3],x[4],x[5]);} };
 struct Vector3r_pickle: bp::pickle_suite{static bp::tuple getinitargs(const Vector3r& x){ return bp::make_tuple(x[0],x[1],x[2]);} };
 struct Vector3i_pickle: bp::pickle_suite{static bp::tuple getinitargs(const Vector3i& x){ return bp::make_tuple(x[0],x[1],x[2]);} };
 struct Vector2r_pickle: bp::pickle_suite{static bp::tuple getinitargs(const Vector2r& x){ return bp::make_tuple(x[0],x[1]);} };
@@ -58,45 +63,24 @@
 
 #undef IDX_CHECK
 
-// automagic converters from sequence (list, tuple, …) to Vector{2,3}{r,i}
-struct custom_Vector3r_from_sequence{
-	custom_Vector3r_from_sequence(){	bp::converter::registry::push_back(&convertible,&construct,bp::type_id<Vector3r>()); }
-	static void* convertible(PyObject* obj_ptr){ if(!PySequence_Check(obj_ptr) || PySequence_Size(obj_ptr)!=3) return 0;	return obj_ptr; }
-	static void construct(PyObject* obj_ptr, bp::converter::rvalue_from_python_stage1_data* data){
-		void* storage=((bp::converter::rvalue_from_python_storage<Vector3r>*)(data))->storage.bytes;
-		new (storage) Vector3r(bp::extract<Real>(PySequence_GetItem(obj_ptr,0)),bp::extract<Real>(PySequence_GetItem(obj_ptr,1)),bp::extract<Real>(PySequence_GetItem(obj_ptr,2)));
-		data->convertible=storage;
-	}
-};
-struct custom_Vector3i_from_sequence{
-	custom_Vector3i_from_sequence(){	bp::converter::registry::push_back(&convertible,&construct,bp::type_id<Vector3i>()); }
-	static void* convertible(PyObject* obj_ptr){ if(!PySequence_Check(obj_ptr) || PySequence_Size(obj_ptr)!=3) return 0;	return obj_ptr; }
-	static void construct(PyObject* obj_ptr, bp::converter::rvalue_from_python_stage1_data* data){
-		void* storage=((bp::converter::rvalue_from_python_storage<Vector3i>*)(data))->storage.bytes;
-		new (storage) Vector3i(bp::extract<int>(PySequence_GetItem(obj_ptr,0)),bp::extract<int>(PySequence_GetItem(obj_ptr,1)),bp::extract<int>(PySequence_GetItem(obj_ptr,2)));
-		data->convertible=storage;
-	}
-};
-struct custom_Vector2r_from_sequence{
-	custom_Vector2r_from_sequence(){	bp::converter::registry::push_back(&convertible,&construct,bp::type_id<Vector2r>()); }
-	static void* convertible(PyObject* obj_ptr){ if(!PySequence_Check(obj_ptr) || PySequence_Size(obj_ptr)!=2) return 0;	return obj_ptr; }
-	static void construct(PyObject* obj_ptr, bp::converter::rvalue_from_python_stage1_data* data){
-		void* storage=((bp::converter::rvalue_from_python_storage<Vector2r>*)(data))->storage.bytes;
-		new (storage) Vector2r(bp::extract<Real>(PySequence_GetItem(obj_ptr,0)),bp::extract<Real>(PySequence_GetItem(obj_ptr,1)));
-		data->convertible=storage;
-	}
-};
-struct custom_Vector2i_from_sequence{
-	custom_Vector2i_from_sequence(){	bp::converter::registry::push_back(&convertible,&construct,bp::type_id<Vector2i>()); }
-	static void* convertible(PyObject* obj_ptr){ if(!PySequence_Check(obj_ptr) || PySequence_Size(obj_ptr)!=2) return 0;	return obj_ptr; }
-	static void construct(PyObject* obj_ptr, bp::converter::rvalue_from_python_stage1_data* data){
-		void* storage=((bp::converter::rvalue_from_python_storage<Vector2i>*)(data))->storage.bytes;
-		new (storage) Vector2i(bp::extract<int>(PySequence_GetItem(obj_ptr,0)),bp::extract<int>(PySequence_GetItem(obj_ptr,1)));
+/* template to define custom converter from sequence/list or approriate length and type, to eigen's Vector
+   - length is stored in VT::RowsAtCompileTime
+	- type is VT::Scalar
+*/
+template<class VT>
+struct custom_VectorAnyAny_from_sequence{
+	custom_VectorAnyAny_from_sequence(){ bp::converter::registry::push_back(&convertible,&construct,bp::type_id<VT>()); }
+	static void* convertible(PyObject* obj_ptr){ if(!PySequence_Check(obj_ptr) || PySequence_Size(obj_ptr)!=VT::RowsAtCompileTime) return 0; return obj_ptr; }
+	static void construct(PyObject* obj_ptr, bp::converter::rvalue_from_python_stage1_data* data){
+		void* storage=((bp::converter::rvalue_from_python_storage<VT>*)(data))->storage.bytes;
+		new (storage) VT;
+		for(size_t i=0; i<VT::RowsAtCompileTime; i++) (*((VT*)storage))[i]=bp::extract<typename VT::Scalar>(PySequence_GetItem(obj_ptr,i));
 		data->convertible=storage;
 	}
 };
 
 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 Vector6r* Vector6r_fromElements(Real v0, Real v1, Real v2, Real v3, Real v4, Real v5){ Vector6r* v(new Vector6r); (*v)<<v0,v1,v2,v3,v4,v5; return v; }
 static Vector3r Matrix3r_diagonal(const Matrix3r& m){ return Vector3r(m.diagonal()); }
 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; }
@@ -125,6 +109,7 @@
 
 EIG_WRAP_METH0(Matrix3r,Zero);
 EIG_WRAP_METH0(Matrix3r,Identity);
+EIG_WRAP_METH0(Vector6r,Zero); EIG_WRAP_METH0(Vector6r,Ones);
 EIG_WRAP_METH0(Vector3r,Zero); EIG_WRAP_METH0(Vector3r,UnitX); EIG_WRAP_METH0(Vector3r,UnitY); EIG_WRAP_METH0(Vector3r,UnitZ); EIG_WRAP_METH0(Vector3r,Ones);
 EIG_WRAP_METH0(Vector3i,Zero); EIG_WRAP_METH0(Vector3i,UnitX); EIG_WRAP_METH0(Vector3i,UnitY); EIG_WRAP_METH0(Vector3i,UnitZ); EIG_WRAP_METH0(Vector3i,Ones);
 EIG_WRAP_METH0(Vector2r,Zero); EIG_WRAP_METH0(Vector2r,UnitX); EIG_WRAP_METH0(Vector2r,UnitY); EIG_WRAP_METH0(Vector2r,Ones);
@@ -146,6 +131,13 @@
 EIG_OP2(Matrix3r,__div__,/,Real) EIG_OP2_INPLACE(Matrix3r,__idiv__,/=,Real)
 EIG_OP2(Matrix3r,__div__,/,int) EIG_OP2_INPLACE(Matrix3r,__idiv__,/=,int)
 
+EIG_OP1(Vector6r,__neg__,-);
+EIG_OP2(Vector6r,__add__,+,Vector6r); EIG_OP2_INPLACE(Vector6r,__iadd__,+=,Vector6r)
+EIG_OP2(Vector6r,__sub__,-,Vector6r); EIG_OP2_INPLACE(Vector6r,__isub__,-=,Vector6r)
+EIG_OP2(Vector6r,__mul__,*,Real) EIG_OP2(Vector6r,__rmul__,*,Real) EIG_OP2_INPLACE(Vector6r,__imul__,*=,Real) EIG_OP2(Vector6r,__div__,/,Real) EIG_OP2_INPLACE(Vector6r,__idiv__,/=,Real)
+EIG_OP2(
+Vector6r,__mul__,*,int) EIG_OP2(Vector6r,__rmul__,*,int) EIG_OP2_INPLACE(Vector6r,__imul__,*=,int) EIG_OP2(Vector6r,__div__,/,int) EIG_OP2_INPLACE(Vector6r,__idiv__,/=,int)
+
 EIG_OP1(Vector3r,__neg__,-);
 EIG_OP2(Vector3r,__add__,+,Vector3r); EIG_OP2_INPLACE(Vector3r,__iadd__,+=,Vector3r)
 EIG_OP2(Vector3r,__sub__,-,Vector3r); EIG_OP2_INPLACE(Vector3r,__isub__,-=,Vector3r)
@@ -175,10 +167,12 @@
 
 	YADE_SET_DOCSTRING_OPTS;
 
-	custom_Vector3r_from_sequence();
-	custom_Vector3i_from_sequence();
-	custom_Vector2r_from_sequence();
-	custom_Vector2i_from_sequence();
+
+	custom_VectorAnyAny_from_sequence<Vector6r>();
+	custom_VectorAnyAny_from_sequence<Vector3r>();
+	custom_VectorAnyAny_from_sequence<Vector3i>();
+	custom_VectorAnyAny_from_sequence<Vector2r>();
+	custom_VectorAnyAny_from_sequence<Vector2i>();
 
 	bp::class_<Matrix3r>("Matrix3","3x3 float matrix.\n\nSupported operations (``m`` is a Matrix3, ``f`` if a float/int, ``v`` is a Vector3): ``-m``, ``m+m``, ``m+=m``, ``m-m``, ``m-=m``, ``m*f``, ``f*m``, ``m*=f``, ``m/f``, ``m/=f``, ``m*m``, ``m*=m``, ``m*v``, ``v*m``, ``m==m``, ``m!=m``.",bp::init<>())
 		.def(bp::init<Matrix3r const &>((bp::arg("m"))))
@@ -239,7 +233,32 @@
 		.def("__setitem__",&Quaternionr_set_item).def("__getitem__",&Quaternionr_get_item)
 		.def("__str__",&Quaternionr_str).def("__repr__",&Quaternionr_str)
 	;
-	bp::class_<Vector3r>("Vector3","3-dimensional float vector.\n\nSupported operatrions (``f`` if a float/int, ``v`` is a Vector3): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*f``, ``f*v``, ``v*=f``, ``v/f``, ``v/=f``, ``v==v``, ``v!=v``, plus operations with ``Matrix3`` and ``Quaternion``.\n\nImplicit conversion from sequence (list,tuple, …) of 3 floats.",bp::init<>())
+	bp::class_<Vector6r>("Vector6","6-dimensional float vector.\n\nSupported operations (``f`` if a float/int, ``v`` is a Vector6): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*f``, ``f*v``, ``v*=f``, ``v/f``, ``v/=f``, ``v==v``, ``v!=v``.\n\nImplicit conversion from sequence (list,tuple, …) of 6 floats.",bp::init<>())
+		.def(bp::init<Vector6r>((bp::arg("other"))))
+		.def("__init__",bp::make_constructor(&Vector6r_fromElements,bp::default_call_policies(),(bp::arg("v0"),bp::arg("v1"),bp::arg("v2"),bp::arg("v3"),bp::arg("v4"),bp::arg("v5"))))
+		.def_pickle(Vector6r_pickle())
+		// properties
+		.add_static_property("Ones",&Vector6r_Ones).add_static_property("Zero",&Vector6r_Zero)
+		//.add_static_property("UnitX",&Vector6r_UnitX).add_static_property("UnitY",&Vector6r_UnitY).add_static_property("UnitZ",&Vector6r_UnitZ)
+		// methods
+		//.def("dot",&Vector6r_dot).def("cross",&Vector6r_cross)
+		.def("norm",&Vector6r::norm).def("squaredNorm",&Vector6r::squaredNorm).def("normalize",&Vector6r::normalize).def("normalized",&Vector6r::normalized)
+		// operators
+		.def("__neg__",&Vector6r__neg__) // -v
+		.def("__add__",&Vector6r__add__Vector6r).def("__iadd__",&Vector6r__iadd__Vector6r) // +, +=
+		.def("__sub__",&Vector6r__sub__Vector6r).def("__isub__",&Vector6r__isub__Vector6r) // -, -=
+		.def("__mul__",&Vector6r__mul__Real).def("__rmul__",&Vector6r__rmul__Real) // f*v, v*f
+		.def("__div__",&Vector6r__div__Real).def("__idiv__",&Vector6r__idiv__Real) // v/f, v/=f
+		.def("__mul__",&Vector6r__mul__int).def("__rmul__",&Vector6r__rmul__int) // f*v, v*f
+		.def("__div__",&Vector6r__div__int).def("__idiv__",&Vector6r__idiv__int) // v/f, v/=f
+		.def(bp::self != bp::self).def(bp::self == bp::self)
+		// specials
+		.def("__len__",&::Vector6r_len).staticmethod("__len__")
+		.def("__setitem__",&::Vector6r_set_item).def("__getitem__",&::Vector6r_get_item)
+		.def("__str__",&::Vector6r_str).def("__repr__",&::Vector6r_str)
+	;
+
+	bp::class_<Vector3r>("Vector3","3-dimensional float vector.\n\nSupported operations (``f`` if a float/int, ``v`` is a Vector3): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*f``, ``f*v``, ``v*=f``, ``v/f``, ``v/=f``, ``v==v``, ``v!=v``, plus operations with ``Matrix3`` and ``Quaternion``.\n\nImplicit conversion from sequence (list,tuple, …) of 3 floats.",bp::init<>())
 		.def(bp::init<Vector3r>((bp::arg("other"))))
 		.def(bp::init<Real,Real,Real>((bp::arg("x"),bp::arg("y"),bp::arg("z"))))
 		.def_pickle(Vector3r_pickle())
@@ -263,7 +282,7 @@
 		.def("__setitem__",&::Vector3r_set_item).def("__getitem__",&::Vector3r_get_item)
 		.def("__str__",&::Vector3r_str).def("__repr__",&::Vector3r_str)
 	;	
-	bp::class_<Vector3i>("Vector3i","3-dimensional integer vector.\n\nSupported operatrions (``i`` if an int, ``v`` is a Vector3i): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*i``, ``i*v``, ``v*=i``, ``v==v``, ``v!=v``.\n\nImplicit conversion from sequence  (list,tuple, …) of 3 integers.",bp::init<>())
+	bp::class_<Vector3i>("Vector3i","3-dimensional integer vector.\n\nSupported operations (``i`` if an int, ``v`` is a Vector3i): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*i``, ``i*v``, ``v*=i``, ``v==v``, ``v!=v``.\n\nImplicit conversion from sequence  (list,tuple, …) of 3 integers.",bp::init<>())
 		.def(bp::init<Vector3i>((bp::arg("other"))))
 		.def(bp::init<int,int,int>((bp::arg("x"),bp::arg("y"),bp::arg("z"))))
 		.def_pickle(Vector3i_pickle())
@@ -284,7 +303,7 @@
 		.def("__setitem__",&::Vector3i_set_item).def("__getitem__",&::Vector3i_get_item)
 		.def("__str__",&::Vector3i_str).def("__repr__",&::Vector3i_str)
 	;	
-	bp::class_<Vector2r>("Vector2","3-dimensional float vector.\n\nSupported operatrions (``f`` if a float/int, ``v`` is a Vector3): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*f``, ``f*v``, ``v*=f``, ``v/f``, ``v/=f``, ``v==v``, ``v!=v``.\n\nImplicit conversion from sequence (list,tuple, …) of 2 floats.",bp::init<>())
+	bp::class_<Vector2r>("Vector2","3-dimensional float vector.\n\nSupported operations (``f`` if a float/int, ``v`` is a Vector3): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*f``, ``f*v``, ``v*=f``, ``v/f``, ``v/=f``, ``v==v``, ``v!=v``.\n\nImplicit conversion from sequence (list,tuple, …) of 2 floats.",bp::init<>())
 		.def(bp::init<Vector2r>((bp::arg("other"))))
 		.def(bp::init<Real,Real>((bp::arg("x"),bp::arg("y"))))
 		.def_pickle(Vector2r_pickle())
@@ -308,7 +327,7 @@
 		.def("__setitem__",&::Vector2r_set_item).def("__getitem__",&::Vector2r_get_item)
 		.def("__str__",&::Vector2r_str).def("__repr__",&::Vector2r_str)
 	;	
-	bp::class_<Vector2i>("Vector2i","2-dimensional integer vector.\n\nSupported operatrions (``i`` if an int, ``v`` is a Vector2i): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*i``, ``i*v``, ``v*=i``, ``v==v``, ``v!=v``.\n\nImplicit conversion from sequence (list,tuple, …) of 2 integers.",bp::init<>())
+	bp::class_<Vector2i>("Vector2i","2-dimensional integer vector.\n\nSupported operations (``i`` if an int, ``v`` is a Vector2i): ``-v``, ``v+v``, ``v+=v``, ``v-v``, ``v-=v``, ``v*i``, ``i*v``, ``v*=i``, ``v==v``, ``v!=v``.\n\nImplicit conversion from sequence (list,tuple, …) of 2 integers.",bp::init<>())
 		.def(bp::init<Vector2i>((bp::arg("other"))))
 		.def(bp::init<int,int>((bp::arg("x"),bp::arg("y"))))
 		.def_pickle(Vector2i_pickle())

=== modified file 'py/tests/wrapper.py'
--- py/tests/wrapper.py	2010-09-30 18:00:41 +0000
+++ py/tests/wrapper.py	2010-10-15 16:49:26 +0000
@@ -32,16 +32,16 @@
 	def testDispatcherCtor(self):
 		# dispatchers take list of their functors in the ctor
 		# same functors are collapsed in one
-		cld1=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_Basic(),Law2_Dem3DofGeom_FrictPhys_Basic()]); self.assert_(len(cld1.functors)==1)
+		cld1=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_CundallStrack(),Law2_Dem3DofGeom_FrictPhys_CundallStrack()]); self.assert_(len(cld1.functors)==1)
 		# two different make two different, right?
-		cld2=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_Basic(),Law2_Dem3DofGeom_CpmPhys_Cpm()]); self.assert_(len(cld2.functors)==2)
-	def testInteractionDispatchersCtor(self):
+		cld2=LawDispatcher([Law2_Dem3DofGeom_FrictPhys_CundallStrack(),Law2_Dem3DofGeom_CpmPhys_Cpm()]); self.assert_(len(cld2.functors)==2)
+	def testInteractionLoopCtor(self):
 		# InteractionLoop takes 3 lists
-		id=InteractionLoop([Ig2_Facet_Sphere_Dem3DofGeom(),Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_Basic()],)
+		id=InteractionLoop([Ig2_Facet_Sphere_Dem3DofGeom(),Ig2_Sphere_Sphere_Dem3DofGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_Dem3DofGeom_FrictPhys_CundallStrack()],)
 		self.assert_(len(id.geomDispatcher.functors)==2)
 		self.assert_(id.geomDispatcher.__class__==IGeomDispatcher)
 		self.assert_(id.physDispatcher.functors[0].__class__==Ip2_FrictMat_FrictMat_FrictPhys)
-		self.assert_(id.lawDispatcher.functors[0].__class__==Law2_Dem3DofGeom_FrictPhys_Basic)
+		self.assert_(id.lawDispatcher.functors[0].__class__==Law2_Dem3DofGeom_FrictPhys_CundallStrack)
 	def testParallelEngineCtor(self):
 		pe=ParallelEngine([InsertionSortCollider(),[BoundDispatcher(),ForceResetter()]])
 		self.assert_(pe.slaves[0].__class__==InsertionSortCollider)