yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2529
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Wed 2010-11-03 17:18:49 +0100
message:
  - Remove scene pointer from TW function parameters
  - sphere-cylinder Ig2 functor (some old uncommited code - still bugged)
modified:
  pkg/common/Cylinder.cpp
  pkg/common/Cylinder.hpp
  pkg/dem/TesselationWrapper.cpp
  pkg/dem/TesselationWrapper.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/common/Cylinder.cpp'
--- pkg/common/Cylinder.cpp	2010-10-20 11:31:49 +0000
+++ pkg/common/Cylinder.cpp	2010-11-03 16:18:49 +0000
@@ -38,30 +38,88 @@
 							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->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
-		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;
-		scm->precompute(state1,state2,scene,c,normal,isNew,shift2,true);
-		return true;
+// 	cerr<<"Ig2_Sphere_ChainedCylinder_CylScGeom::go"<<endl;
+	const State* sphereSt=YADE_CAST<const State*>(&state1);
+	const ChainedState* cylinderSt=YADE_CAST<const ChainedState*>(&state2);
+	ChainedCylinder *cylinder=YADE_CAST<ChainedCylinder*>(cm2.get());
+	Sphere *sphere=YADE_CAST<Sphere*>(cm1.get());
+	assert(sphereSt && cylinderSt && cylinder && sphere);
+	if (cylinderSt->chains[cylinderSt->chainNumber].size()==(cylinderSt->rank+1)) {cerr << "last cylinder - ignored"<<endl; return false;}
+
+// 	int i1=cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1];
+// 	cerr << "i1 : "<<i1<< " i2: "<<cylinderSt->rank<<" pos: "<<Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1],scene)->state->pos<< " - "<<cylinderSt->pos<<endl;
+
+	//FIXME : definition of segment in next line breaks periodicity
+	Vector3r segment = Body::byId(cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1],scene)->state->pos-cylinderSt->pos;
+	Vector3r branch = sphereSt->pos-cylinderSt->pos-shift2;
+	bool isNew = !c->geom;
+
+	//Check position of projection on cylinder axis
+// 	if ((segment.dot(branch)>(segment.dot(segment)/*+interactionDetectionFactor*cylinder->radius*/) && !c->isReal())) return (false);//position _after_ end of cylinder
+	if ((segment.dot(branch)>(segment.dot(segment)/*+interactionDetectionFactor*cylinder->radius*/) && isNew)) return (false);//position _after_ end of cylinder
+// 	if ((segment.dot(branch)>(segment.dot(segment)/*+interactionDetectionFactor*cylinder->radius*/) && c->isReal())) cerr<<"pb1"<<endl;
+
+	Real length = segment.norm();
+	Vector3r direction = segment/length;
+	Real dist = direction.dot(branch);
+// 	if ((dist<-interactionDetectionFactor*cylinder->radius) && !c->isReal()) return (false);//position _before_ start of cylinder
+	if ((dist<-interactionDetectionFactor*cylinder->radius) && isNew) return (false);//position _before_ start of cylinder
+		
+	//Check sphere-cylinder distance
+	Vector3r projectedP = cylinderSt->pos+shift2 + direction*dist;
+	branch = projectedP-sphereSt->pos;
+// 	if ((branch.squaredNorm()>(pow(sphere->radius+cylinder->radius,2))) && !c->isReal()) return (false);
+	if ((branch.squaredNorm()>(pow(sphere->radius+cylinder->radius,2))) && isNew) return (false);
+
+	shared_ptr<CylScGeom> scm;
+	if(!isNew) scm=YADE_PTR_CAST<CylScGeom>(c->geom);
+	else { scm=shared_ptr<CylScGeom>(new CylScGeom()); c->geom=scm; }
+	
+	scm->radius1=sphere->radius;
+	scm->radius2=cylinder->radius;
+	scm->id3=cylinderSt->chains[cylinderSt->chainNumber][cylinderSt->rank+1];
+	scm->start=cylinderSt->pos+shift2; scm->end=scm->start+segment;
+
+	//FIXME : there should be other checks without distanceFactor?
+	if (dist<0.0) {//We have sphere-node contact
+		Vector3r normal=(cylinderSt->pos+shift2)-sphereSt->pos;
+		Real norm=normal.norm();
+		scm->onNode=true; scm->relPos=0;
+		scm->penetrationDepth=sphere->radius+cylinder->radius-norm;
+		scm->contactPoint=sphereSt->pos+(sphere->radius-0.5*scm->penetrationDepth)*normal;
+		scm->precompute(state1,state2,scene,c,normal/norm,isNew,shift2,true);//use sphere-sphere precompute (a node is a sphere)
+	} else {//we have sphere-cylinder contact
+		scm->onNode=false; scm->relPos=dist/length;
+		Real norm=branch.norm();
+		Vector3r normal=branch/norm;
+		scm->penetrationDepth= sphere->radius+cylinder->radius-norm;
+		if (dist>length) {
+			scm->penetrationDepth=sphere->radius+cylinder->radius-(cylinderSt->pos+segment-sphereSt->pos).norm();
+			//FIXME : handle contact jump on next element
+		}
+		scm->contactPoint = sphereSt->pos+scm->normal*(sphere->radius-0.5*scm->penetrationDepth);
+
+		//FIXME : replace the block below with smart use of shift2=shift2 + cylinder_spin.cross(branch) - doesn't compile currently
+
+		//precompute stuff manually (sphere-sphere precompute doesn't apply here if we want to avoid ratcheting
+		/*
+		if(!isNew) {
+			scm->orthonormal_axis = scm->normal.cross(normal);
+			Real angle = scene->dt*0.5*scm->normal.dot(sphereSt->angVel + cylinderSt->angVel);
+			scm->twist_axis = angle*normal;}
+		else scm->twist_axis=scm->orthonormal_axis=Vector3r::Zero();
+		//Update contact normal
+		scm->normal=normal;
+		Vector3r c1x =  sphere->radius*normal;
+		Vector3r c2x = -cylinder->radius*normal;
+		Vector3r relativeVelocity = (cylinderSt->vel+cylinderSt->angVel.cross(c2x)) - (sphereSt->vel+sphereSt->angVel.cross(c1x));
+		//keep the shear part only
+		relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
+		scm->shearInc = relativeVelocity*scene->dt;
+		*/
+
 	}
-	return false;
+	return true;
 }
 
 
@@ -73,7 +131,9 @@
 								const bool& force,
 								const shared_ptr<Interaction>& c)
 {
-	return go(cm1,cm2,state2,state1,-shift2,force,c);
+ 	cerr<<"Ig2_Sphere_ChainedCylinder_CylScGeom::goReverse"<<endl;
+	c->swapOrder();
+	return go(cm2,cm1,state2,state1,-shift2,force,c);
 }
 
 
@@ -95,7 +155,7 @@
 	if (!pChain1 || !pChain2) {
 		cerr <<"cast failed8567"<<endl;
 	}
-	const bool revert = (pChain2->rank-pChain1->rank == -1);
+	const bool revert = ((int) pChain2->rank- (int) pChain1->rank == -1);
 	const ChainedState& bchain1 = revert? *pChain2 : *YADE_CAST<const ChainedState*>(&state1);
 	const ChainedState& bchain2 = revert? *pChain1 : *pChain2;
 	if (bchain2.rank-bchain1.rank != 1) {/*cerr<<"Mutual contacts in same chain between not adjacent elements, not handled*/ return false;}

=== modified file 'pkg/common/Cylinder.hpp'
--- pkg/common/Cylinder.hpp	2010-10-13 16:23:08 +0000
+++ pkg/common/Cylinder.hpp	2010-11-03 16:18:49 +0000
@@ -36,7 +36,7 @@
 
 	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"))
-  		((Quaternionr,chainedOrientation,Quaternionr::Identity(),,"Orientation of node-to-node vector"))
+  		((Quaternionr,chainedOrientation,Quaternionr::Identity(),,"Deviation of node1 orientation from node-to-node vector"))
 		,createIndex();/*ctor*/
 // 		state=shared_ptr<ChainedState>(new ChainedState);
 
@@ -48,7 +48,11 @@
 	public:
 		virtual ~CylScGeom ();
 	YADE_CLASS_BASE_DOC_ATTRS_CTOR(CylScGeom,ScGeom,"Geometry of a cylinder-sphere contact.",
-		((bool,onNode,false,,"contact on node")),
+		((bool,onNode,false,,"contact on node?"))
+		((Vector3r,start,Vector3r::Zero(),,"position of 1st node |yupdate|"))
+		((Vector3r,end,Vector3r::Zero(),,"position of 2nd node |yupdate|"))
+		((Body::id_t,id3,0,,"id of next chained cylinder |yupdate|"))
+		((Real,relPos,0,,"position of the contact on the cylinder (0: node-, 1:node+) |yupdate|")),
 		createIndex(); /*ctor*/
 	);
 	REGISTER_CLASS_INDEX(CylScGeom,ScGeom);
@@ -57,21 +61,21 @@
 
 class ChainedState: public State{
 	public:
-		static vector<vector<int> > chains;
+		static vector<vector<Body::id_t> > chains;
 		static unsigned int currentChain;
-		vector<int> barContacts;
-		vector<int> nodeContacts;
+		vector<Body::id_t> barContacts;
+		vector<Body::id_t> nodeContacts;
 
 		virtual ~ChainedState ();
-		void addToChain(int bodyId) {
+		void addToChain(Body::id_t bodyId) {
 			if (chains.size()<=currentChain) chains.resize(currentChain+1);
 			chainNumber=currentChain;
  			rank=chains[currentChain].size();
- 			chains[currentChain].push_back(rank);}
+ 			chains[currentChain].push_back(bodyId);}
 
 	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"))
+ 		((unsigned int,rank,0,,"rank in the chain"))
+ 		((unsigned int,chainNumber,0,,"chain id"))
 		,
 		/* additional initializers */
 /*			((pos,se3.position))

=== modified file 'pkg/dem/TesselationWrapper.cpp'
--- pkg/dem/TesselationWrapper.cpp	2010-11-02 16:23:44 +0000
+++ pkg/dem/TesselationWrapper.cpp	2010-11-03 16:18:49 +0000
@@ -11,6 +11,7 @@
 #include<yade/extra/boost_python_len.hpp>
 #include<yade/pkg-dem/Shop.hpp>
 #include"TesselationWrapper.hpp"
+#include <yade/lib-triangulation/Timer.h>
 
 YADE_PLUGIN((TesselationWrapper));
 YADE_REQUIRE_FEATURE(CGAL)
@@ -154,10 +155,14 @@
 //  facet_it = Tes->Triangulation().finite_edges_end ();
 }
 
-void TesselationWrapper::insertSceneSpheres(const Scene* scene)
+void TesselationWrapper::insertSceneSpheres()
 {
-	const shared_ptr<BodyContainer>& bodies = scene->bodies;
+	Scene* scene=Omega::instance().getScene().get();
+// 	Real_timer clock;
+//         clock.start();
+        const shared_ptr<BodyContainer>& bodies = scene->bodies;
 	build_triangulation_with_ids(bodies, *this);
+// 	clock.top("Triangulation");
 }
 
 double TesselationWrapper::Volume(unsigned int id) {return Tes->Volume(id);}

=== modified file 'pkg/dem/TesselationWrapper.hpp'
--- pkg/dem/TesselationWrapper.hpp	2010-11-02 17:40:46 +0000
+++ pkg/dem/TesselationWrapper.hpp	2010-11-03 16:18:49 +0000
@@ -27,7 +27,7 @@
  *O.load("test.xml")
  *O.run(100) //for unknown reasons, this procedure crashes at iteration 0
  *TW=TesselationWrapper()
- *TW.triangulate(O._sceneObj()) //compute regular Delaunay triangulation, don't construct tesselation
+ *TW.triangulate() //compute regular Delaunay triangulation, don't construct tesselation
  *TW.computeVolumes() //will silently tesselate the packing
  *TW.volume(10) //get volume associated to sphere of id 10
  * 
@@ -49,7 +49,7 @@
     	/// Insert a sphere, "id" will be used by some getters to retrieve spheres
     	bool insert(double x, double y, double z, double rad, unsigned int id);
 	/// A faster version of insert, inserting all spheres in scene
-	void insertSceneSpheres(const Scene* scene);
+	void insertSceneSpheres();
 	/// Move one sphere to the new position (x,y,z) and maintain triangulation (invalidates the tesselation)
 	bool move (double x, double y, double z, double rad, unsigned int id);
 	
@@ -93,7 +93,7 @@
 	CGT::Finite_edges_iterator facet_it;
 	MicroMacroAnalyser mma;
 	
-	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(TesselationWrapper,GlobalEngine,"Handle the triangulation of spheres in a scene, build tesselation on request, and give access to computed quantities : currently volume and porosity of each Voronoï sphere. Example script :\n tt=TriaxialTest()\ntt.generate('test.xml')\nO.load('test.xml')\nO.run(100) //for unknown reasons, this procedure crashes at iteration 0\nTW=TesselationWrapper()\nTW.triangulate(O._sceneObj()) //compute regular Delaunay triangulation, don't construct tesselation\nTW.computeVolumes() //will silently tesselate the packing\nTW.volume(10) //get volume associated to sphere of id 10",
+	YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(TesselationWrapper,GlobalEngine,"Handle the triangulation of spheres in a scene, build tesselation on request, and give access to computed quantities : currently volume and porosity of each Voronoï sphere. Example script :\n tt=TriaxialTest()\ntt.generate('test.xml')\nO.load('test.xml')\nO.run(100) //for unknown reasons, this procedure crashes at iteration 0\nTW=TesselationWrapper()\nTW.triangulate() //compute regular Delaunay triangulation, don't construct tesselation\nTW.computeVolumes() //will silently tesselate the packing\nTW.volume(10) //get volume associated to sphere of id 10",
 	((unsigned int,n_spheres,0,,"|ycomp|"))
 	,/*ctor*/
   	Tes = new CGT::Tesselation;