yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2517
committer: bchareyre <bchareyre@dt-rv020>
branch nick: yade
timestamp: Wed 2010-10-27 12:14:05 +0200
message:
  - Convert ScGeom to ScGeom6D in Ig2_Sphere_Sphere_ScGeom6D::go(), using ScGeom::operator=
  - Optimize distance check for existing interactions in Ig2_Sphere_Sphere_ScGeom::go()
  - Implement Ig2_Box_Sphere_ScGeom6D (same conversion to ScGeom6D)
  - Adapt CohesiveTriaxialTest.
modified:
  pkg/dem/CohesiveTriaxialTest.cpp
  pkg/dem/Ig2_Box_Sphere_ScGeom.cpp
  pkg/dem/Ig2_Box_Sphere_ScGeom.hpp
  pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp
  pkg/dem/ScGeom.cpp
  pkg/dem/ScGeom.hpp


--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk

Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription

=== modified file 'pkg/dem/CohesiveTriaxialTest.cpp'
--- pkg/dem/CohesiveTriaxialTest.cpp	2010-10-25 14:52:21 +0000
+++ pkg/dem/CohesiveTriaxialTest.cpp	2010-10-27 10:14:05 +0000
@@ -226,6 +226,7 @@
 	physics->frictionAngle		= sphereFrictionDeg * Mathr::PI/180.0;
 	physics->shearCohesion = shearCohesion;
 	physics->normalCohesion = normalCohesion;
+	physics->momentRotationLaw = 1;
 	
 	if((!dynamic) && (!boxWalls))
 	{
@@ -274,6 +275,7 @@
 	physics->frictionAngle		= boxFrictionDeg * Mathr::PI/180.0;
 	physics->shearCohesion = 0;
 	physics->normalCohesion = 0;
+	physics->momentRotationLaw = 0;
 	
 	aabb->color		= Vector3r(1,0,0);
 	
@@ -293,7 +295,7 @@
 	shared_ptr<IGeomDispatcher> interactionGeometryDispatcher(new IGeomDispatcher);
 	shared_ptr<IGeomFunctor> s1(new Ig2_Sphere_Sphere_ScGeom6D);
 	interactionGeometryDispatcher->add(s1);
-	shared_ptr<IGeomFunctor> s2(new Ig2_Box_Sphere_ScGeom);
+	shared_ptr<IGeomFunctor> s2(new Ig2_Box_Sphere_ScGeom6D);
 	interactionGeometryDispatcher->add(s2);
 
 	shared_ptr<Ip2_2xCohFrictMat_CohFrictPhys> cohesiveFrictionalRelationships = shared_ptr<Ip2_2xCohFrictMat_CohFrictPhys> (new Ip2_2xCohFrictMat_CohFrictPhys);

=== modified file 'pkg/dem/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Ig2_Box_Sphere_ScGeom.cpp	2010-10-16 18:31:17 +0000
+++ pkg/dem/Ig2_Box_Sphere_ScGeom.cpp	2010-10-27 10:14:05 +0000
@@ -147,12 +147,14 @@
 		scm->radius1 = s->radius;
 		scm->radius2 = s->radius;
 		c->geom = scm;
-		scm->precompute(state1,state2,scene,c,-cOnBox_sphere,isNew,shift2,true);
+		//FIXME : NOT TESTED : do we precompute correctly if boxes are moving?
+		scm->precompute(state1,state2,scene,c,-cOnBox_sphere,isNew,shift2,false);
 	}
 	return true;
 }
 
 
+
 bool Ig2_Box_Sphere_ScGeom::goReverse(	const shared_ptr<Shape>& cm1,
 						const shared_ptr<Shape>& cm2,
 						const State& state1,
@@ -166,3 +168,45 @@
 }
 
 YADE_PLUGIN((Ig2_Box_Sphere_ScGeom));
+
+
+#ifdef YADE_DEVIRT_FUNCTORS
+bool Ig2_Box_Sphere_ScGeom6D::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_Box_Sphere_ScGeom6D::go, use getStaticFunctorPtr and call that function instead."); }
+bool Ig2_Box_Sphere_ScGeom::goStatic(IGeomFunctor* 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){
+#else
+bool Ig2_Box_Sphere_ScGeom6D::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
+	bool isNew = !c->geom;
+	if (Ig2_Box_Sphere_ScGeom::go(cm1,cm2,state1,state2,shift2,force,c)) {
+		if (isNew) {//generate a 6DOF interaction from the 3DOF one generated by Ig2_Box_Sphere_ScGeom
+			shared_ptr<ScGeom6D> sc (new ScGeom6D());
+			*(YADE_PTR_CAST<ScGeom>(sc)) = *(YADE_PTR_CAST<ScGeom>(c->geom));
+			c->geom=sc;}
+		YADE_PTR_CAST<ScGeom6D>(c->geom)->precomputeRotations(state1,state2,isNew,false);
+		return true;}
+	else return false;
+}
+
+
+
+bool Ig2_Box_Sphere_ScGeom6D::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)
+{
+	c->swapOrder();
+	return go(cm2,cm1,state2,state1,-shift2,force,c);
+}
+
+YADE_PLUGIN((Ig2_Box_Sphere_ScGeom6D));

=== modified file 'pkg/dem/Ig2_Box_Sphere_ScGeom.hpp'
--- pkg/dem/Ig2_Box_Sphere_ScGeom.hpp	2010-10-13 16:23:08 +0000
+++ pkg/dem/Ig2_Box_Sphere_ScGeom.hpp	2010-10-27 10:14:05 +0000
@@ -44,3 +44,33 @@
 };
 REGISTER_SERIALIZABLE(Ig2_Box_Sphere_ScGeom);
 
+class Ig2_Box_Sphere_ScGeom6D : public Ig2_Box_Sphere_ScGeom
+{
+	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_Box_Sphere_ScGeom::goStatic; }
+		static bool goStatic(IGeomFunctor* 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(Ig2_Box_Sphere_ScGeom6D,Ig2_Box_Sphere_ScGeom,"Create an interaction geometry :yref:`ScGeom6D` from :yref:`Box` and :yref:`Sphere`")
+	FUNCTOR2D(Box,Sphere);
+	DEFINE_FUNCTOR_ORDER_2D(Box,Sphere);
+};
+REGISTER_SERIALIZABLE(Ig2_Box_Sphere_ScGeom6D);
+

=== modified file 'pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp'
--- pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp	2010-10-20 17:46:55 +0000
+++ pkg/dem/Ig2_Sphere_Sphere_ScGeom.cpp	2010-10-27 10:14:05 +0000
@@ -22,28 +22,26 @@
 {
 #endif
 	const Se3r& se31=state1.se3; const Se3r& se32=state2.se3;
-	Sphere *s1=static_cast<Sphere*>(cm1.get()), *s2=static_cast<Sphere*>(cm2.get());
+	const 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->geom;
-		if(!isNew) scm=YADE_PTR_CAST<ScGeom>(c->geom);
-		else { scm=shared_ptr<ScGeom>(new ScGeom()); c->geom=scm; }
-
-		Real norm=normal.norm(); normal/=norm; // normal is unit vector now
+	if (!c->isReal() && !force) {//don't fast-check distance if geometry will be updated anyway
+		Real penetrationDepthSq=pow(interactionDetectionFactor*(s1->radius+s2->radius),2) - normal.squaredNorm();
+		if (penetrationDepthSq<0) return false;}
+	shared_ptr<ScGeom> scm;
+	bool isNew = !c->geom;
+	if(!isNew) scm=YADE_PTR_CAST<ScGeom>(c->geom);
+	else { scm=shared_ptr<ScGeom>(new ScGeom()); c->geom=scm; }
+	Real norm=normal.norm(); normal/=norm; // normal is unit vector now
 #ifdef YADE_DEBUG
-		if(norm==0) throw runtime_error(("Zero distance between spheres #"+lexical_cast<string>(c->getId1())+" and #"+lexical_cast<string>(c->getId2())+".").c_str());
+	if(norm==0) throw runtime_error(("Zero distance between spheres #"+lexical_cast<string>(c->getId1())+" and #"+lexical_cast<string>(c->getId2())+".").c_str());
 #endif
-		Real penetrationDepth=s1->radius+s2->radius-norm;
-		scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
-		scm->penetrationDepth=penetrationDepth;
-		scm->radius1=s1->radius;
-		scm->radius2=s2->radius;
-		scm->precompute(state1,state2,scene,c,normal,isNew,shift2,avoidGranularRatcheting);
-		return true;
-	}
-	return false;
+	Real penetrationDepth=s1->radius+s2->radius-norm;
+	scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
+	scm->penetrationDepth=penetrationDepth;
+	scm->radius1=s1->radius;
+	scm->radius2=s2->radius;
+	scm->precompute(state1,state2,scene,c,normal,isNew,shift2,avoidGranularRatcheting);
+	return true;
 }
 
 bool Ig2_Sphere_Sphere_ScGeom::goReverse(	const shared_ptr<Shape>& cm1,
@@ -73,8 +71,11 @@
 #endif
 	bool isNew = !c->geom;
 	if (Ig2_Sphere_Sphere_ScGeom::go(cm1,cm2,state1,state2,shift2,force,c)){//the 3 DOFS from ScGeom are updated here
-		shared_ptr<ScGeom6D> scm=YADE_PTR_CAST<ScGeom6D>(c->geom);
-		if (updateRotations) scm->precomputeRotations(state1,state2,isNew,creep);
+ 		if (isNew) {//generate a 6DOF interaction from the 3DOF one generated by Ig2_Sphere_Sphere_ScGeom
+			shared_ptr<ScGeom6D> sc (new ScGeom6D());
+			*(YADE_PTR_CAST<ScGeom>(sc)) = *(YADE_PTR_CAST<ScGeom>(c->geom));
+			c->geom=sc;}
+		if (updateRotations) YADE_PTR_CAST<ScGeom6D>(c->geom)->precomputeRotations(state1,state2,isNew,creep);
 		return true;
 	}
 	else return false;

=== modified file 'pkg/dem/ScGeom.cpp'
--- pkg/dem/ScGeom.cpp	2010-10-26 16:31:33 +0000
+++ pkg/dem/ScGeom.cpp	2010-10-27 10:14:05 +0000
@@ -95,12 +95,17 @@
 		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.
 // add -DYADE_SCGEOM_DEBUG to CXXFLAGS to enable this piece or just do
-//   #define YADE_SCGEOM_DEBUG (but do not commit with that enabled in the code)
+// #define YADE_SCGEOM_DEBUG //(but do not commit with that enabled in the code)
 #ifdef YADE_SCGEOM_DEBUG
 		if (isnan(aa.angle())) {
 			cerr<<"NaN angle found in angleAxisr(q), for quaternion "<<delta<<", after quaternion product"<<endl;
-			cerr<<"rbp1.ori * (initialOrientation1.conjugate())) * (initialOrientation2 * (rbp2.ori.conjugate()) is:"<<endl;
-			cerr<<rbp1.ori<<" * "<<initialOrientation1.conjugate()<<" * "<<initialOrientation2<<" * "<<rbp2.ori.conjugate()<<endl<<"with sub-products :"<<endl<<rbp1.ori * (initialOrientation1.conjugate())<<" * "<<initialOrientation2 * (rbp2.ori.conjugate())<<endl;
+			cerr<<"rbp1.ori * (initialOrientation1.conjugate())) * (initialOrientation2 * (rbp2.ori.conjugate()) with quaternions :"<<endl;
+			cerr<<rbp1.ori<<" * "<<initialOrientation1<<" * "<<initialOrientation2<<" * "<<rbp2.ori<<endl<<" and sub-products :"<<endl<<rbp1.ori * (initialOrientation1.conjugate())<<" * "<<initialOrientation2 * (rbp2.ori.conjugate())<<endl;
+			cerr <<"q.w (before normalization) "<<delta.w();
+			delta.normalize();
+			cerr <<"q.w (after) "<<delta.w()<<endl;
+			AngleAxisr bb(delta);
+			cerr<<delta<<" "<<bb.angle()<<endl;
 		}
 #else
 		if (isnan(aa.angle())) aa.angle()=0;

=== modified file 'pkg/dem/ScGeom.hpp'
--- pkg/dem/ScGeom.hpp	2010-10-20 15:38:06 +0000
+++ pkg/dem/ScGeom.hpp	2010-10-27 10:14:05 +0000
@@ -25,6 +25,12 @@
 		// inherited from GenericSpheresContact: Vector3r& normal; 
 		Real &radius1, &radius2;
 		virtual ~ScGeom();
+		inline ScGeom& operator= (const ScGeom& source){
+			normal=source.normal; contactPoint=source.contactPoint;
+			twist_axis=source.twist_axis; orthonormal_axis=source.orthonormal_axis;
+			radius1=source.radius1; radius2=source.radius2;
+			penetrationDepth=source.penetrationDepth; shearInc=source.shearInc;
+			return *this;}
 		
 		//!precompute values of shear increment and interaction rotation data. Update contact normal to the vurrentNormal value. Precondition : the value of normal is not updated outside (and before) this function.
 		void precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r& currentNormal, bool isNew, const Vector3r& shift2, bool avoidGranularRatcheting=true);
@@ -61,7 +67,7 @@
 
 		void precomputeRotations(const State& rbp1, const State& rbp2, bool isNew, bool creep=false);
 		void initRotations(const State& rbp1, const State& rbp2);
-	
+  		
 		YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom6D,ScGeom,"Class representing :yref:`geometry<IGeom>` of two :yref:`bodies<Body>` in contact. The contact has 6 DOFs (normal, 2×shear, twist, 2xbending) and uses :yref:`ScGeom` incremental algorithm for updating shear.",
 		((Quaternionr,orientationToContact1,Quaternionr(1.0,0.0,0.0,0.0),,""))
 		((Quaternionr,orientationToContact2,Quaternionr(1.0,0.0,0.0,0.0),,""))