yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-01-20 19:30:38 +0100 (Tue, 20 Jan 2009)
New Revision: 1630

Added:
   trunk/scripts/test-sphere-facet.py
Modified:
   trunk/gui/py/PythonUI.cpp
   trunk/gui/py/PythonUI_rc.py
   trunk/gui/qt3/GLSimulationPlayerViewer.cpp
   trunk/gui/qt3/QtGUI.cpp
   trunk/lib/opengl/GLUtils.hpp
   trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp
   trunk/pkg/common/RenderingEngine/GLDrawGeometricalModel/GLDrawSphere.cpp
   trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
   trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
   trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
   trunk/pkg/dem/SConscript
Log:
1. Add algorithm for pushing sphere always to the same side when in contact with facet, even under extreme stress where the sphere geometrically passes to the other side.
2. Add bool SpheresContactGeometry::initContactOnPositiveFacetSide for that
3. Add test/demo script for that.
4. Re-enabled GLDrawSphere::glutUse (cosurgi?)
5. Renamed yade.runtime.args to argv, to be consistent with sys.argv
6. Allow passing sys.argv (arguments to scripts) through QtGUI as well



Modified: trunk/gui/py/PythonUI.cpp
===================================================================
--- trunk/gui/py/PythonUI.cpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/gui/py/PythonUI.cpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -90,7 +90,7 @@
 			PYTHON_DEFINE_STRING("script",runScript);
 			PYTHON_DEFINE_BOOL("stopAfter",stopAfter);
 			PYTHON_DEFINE_BOOL("nonInteractive",nonInteractive);
-			{ ostringstream oss; oss<<"yade.runtime.args=["; if(scriptArgs.size()>0){ FOREACH(string s, scriptArgs) oss<<"'"<<s<<"',"; } oss<<"]"; PyRun_SimpleString(oss.str().c_str()); }
+			{ ostringstream oss; oss<<"yade.runtime.argv=["; if(scriptArgs.size()>0){ FOREACH(string s, scriptArgs) oss<<"'"<<s<<"',"; } oss<<"]"; PyRun_SimpleString(oss.str().c_str()); }
 		#undef PYTHON_DEFINE_STRING
 		#undef PYTHON_DEFINE_BOOL
 		execScript(PREFIX "/lib/yade" SUFFIX "/gui/PythonUI_rc.py");

Modified: trunk/gui/py/PythonUI_rc.py
===================================================================
--- trunk/gui/py/PythonUI_rc.py	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/gui/py/PythonUI_rc.py	2009-01-20 18:30:38 UTC (rev 1630)
@@ -28,6 +28,8 @@
 srv=yade.PythonTCPServer.PythonTCPServer(minPort=9000)
 yade.runtime.cookie=srv.server.cookie
 sys.stdout.flush()
+runtime.argv=[runtime.script]+runtime.argv
+sys.argv=runtime.argv # could be [] as well
 
 ## run simulation if requested from the command line
 if runtime.simulation:
@@ -57,9 +59,8 @@
 	import time;
 	while True: time.sleep(1)
 else:
-	sys.argv=runtime.args # could be [] as well
+	sys.argv[0]=['<embedded python interpreter>']
 	from IPython.Shell import IPShellEmbed
-	sys.argv=['<embedded python interpreter>']
 	ipshell = IPShellEmbed(banner=r"""__   __    ____          ____                      _      
 \ \ / /_ _|  _ \  ___   / ___|___  _ __  ___  ___ | | ___ 
  \ V / _` | | | |/ _ \ | |   / _ \| '_ \/ __|/ _ \| |/ _ \ 

Modified: trunk/gui/qt3/GLSimulationPlayerViewer.cpp
===================================================================
--- trunk/gui/qt3/GLSimulationPlayerViewer.cpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/gui/qt3/GLSimulationPlayerViewer.cpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -149,8 +149,8 @@
 		}
 		xyzNames.sort();
 	} else { /* load from sqlite database */
+		LOG_INFO("Opening sqlite database `"<<fileName<<"'");
 		con=shared_ptr<sqlite3x::sqlite3_connection>(new sqlite3x::sqlite3_connection(fileName));
-		LOG_DEBUG("Opened sqlite db "<<fileName);
 		if(0==con->executeint("select count(*) from sqlite_master where name='meta';")){ LOG_ERROR("Database doesn't have the 'meta' table."); return; }
 		if(0==con->executeint("select count(*) from sqlite_master where name='records';")){ LOG_ERROR("Database doesn't have the 'records' table."); return; }
 		simPlayer->pushMessage("Database OK, loading simulation...");

Modified: trunk/gui/qt3/QtGUI.cpp
===================================================================
--- trunk/gui/qt3/QtGUI.cpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/gui/qt3/QtGUI.cpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -54,8 +54,15 @@
 		#endif
 		else optind--;
 	}
-	for (int index=optind+1; index<argc; index++) LOG_ERROR("Unprocessed non-option argument: `"<<argv[index]<<"'");
-
+	for (int index=optind+1; index<argc; index++) {
+		#ifdef EMBED_PYTHON
+			LOG_DEBUG("Adding script parameter `"<<argv[index]<<"' from the command line.");
+			PythonUI::scriptArgs.push_back(string(argv[index]));
+			if(!PythonUI::runScript.empty()) LOG_WARN("Got parameter `"<<argv[index]<<"', but no .py script to be run!");
+		#else
+			LOG_ERROR("Unprocessed non-option argument: `"<<argv[index]<<"'");
+		#endif
+	}
 	XInitThreads();
    QApplication app(argc,argv);
 	mainWindow=new YadeQtMainWindow();

Modified: trunk/lib/opengl/GLUtils.hpp
===================================================================
--- trunk/lib/opengl/GLUtils.hpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/lib/opengl/GLUtils.hpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -11,7 +11,6 @@
 #include<string>
 
 struct GLUtils{
-
 	static void GLDrawArrow(const Vector3r& from, const Vector3r& to, const Vector3r& color=Vector3r(1,1,1)){
 		glEnable(GL_LIGHTING); glColor3v(color); qglviewer::Vec a(from[0],from[1],from[2]),b(to[0],to[1],to[2]); QGLViewer::drawArrow(a,b);	
 	}

Modified: trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/pkg/common/Engine/DeusExMachina/GravityEngines.hpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -16,10 +16,9 @@
 		virtual ~GravityEngine(){};
 		virtual void applyCondition(MetaBody*);
 	protected :
-		virtual void registerAttributes(){REGISTER_ATTRIBUTE(gravity);}
 	NEEDS_BEX("Force");
-	REGISTER_CLASS_NAME(GravityEngine);
-	REGISTER_BASE_CLASS_NAME(DeusExMachina);
+	REGISTER_ATTRIBUTES(DeusExMachina,(gravity));
+	REGISTER_CLASS_AND_BASE(GravityEngine,DeusExMachina);
 };
 REGISTER_SERIALIZABLE(GravityEngine);
 

Modified: trunk/pkg/common/RenderingEngine/GLDrawGeometricalModel/GLDrawSphere.cpp
===================================================================
--- trunk/pkg/common/RenderingEngine/GLDrawGeometricalModel/GLDrawSphere.cpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/pkg/common/RenderingEngine/GLDrawGeometricalModel/GLDrawSphere.cpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -37,7 +37,7 @@
 
 void GLDrawSphere::go(const shared_ptr<GeometricalModel>& gm, const shared_ptr<PhysicalParameters>& ph,bool wire)
 {
-	/*if(glutUse){
+	if(glutUse){
 		Real r= (static_cast<Sphere*>(gm.get()))->radius;
 		glMaterialv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, Vector3f(gm->diffuseColor[0],gm->diffuseColor[1],gm->diffuseColor[2]));
 		glColor3v(gm->diffuseColor);
@@ -46,7 +46,7 @@
 			else glutSolidSphere(r,glutSlices,glutStacks);
 		if(glutNormalize) glPopAttrib();
 		return;
-	}*/
+	}
 
 	//first = true;
 	if (first)

Modified: trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp
===================================================================
--- trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/pkg/dem/DataClass/InteractionGeometry/SpheresContactGeometry.hpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -32,9 +32,14 @@
 			contactPoint;
 		Real radius1,radius2,penetrationDepth;
 
+
+		// begin abusive storage
 		bool hasNormalViscosity;
 		Real NormalViscisity;
 		Real NormalRelativeVelocity;
+		//! Whether the original contact was on the positive or negative facet side; this is to permit repulsion to the right side even if the sphere passes, under extreme pressure, geometrically to the other facet's side. This is used only in InteractingFacet2IteractingSphere4SpheresContactGeometry.
+		bool initContactOnPositiveFacetSide;
+		// end abusive storage
 
 		bool hasShear; // whether the exact rotation code is being used -- following fields are needed for that
 		//! positions and orientations of both spheres -- must be updated at every iteration
@@ -93,28 +98,26 @@
 
 		SpheresContactGeometry():contactPoint(Vector3r::ZERO),radius1(0),radius2(0),hasShear(false),pos1(Vector3r::ZERO),pos2(Vector3r::ZERO),ori1(Quaternionr::IDENTITY),ori2(Quaternionr::IDENTITY),cp1rel(Quaternionr::IDENTITY),cp2rel(Quaternionr::IDENTITY),d1(0),d2(0),d0(0),initRelOri12(Quaternionr::IDENTITY){createIndex();}
 		virtual ~SpheresContactGeometry();
-	protected :
-		virtual void registerAttributes(){
-			REGISTER_ATTRIBUTE(normal);
-			REGISTER_ATTRIBUTE(contactPoint);
-			REGISTER_ATTRIBUTE(radius1);
-			REGISTER_ATTRIBUTE(radius2);
-			REGISTER_ATTRIBUTE(contactPoint); // to allow access from python
+	REGISTER_ATTRIBUTES(/* no attributes from base class */,
+			(normal)
+			(contactPoint)
+			(radius1)
+			(radius2)
+			(contactPoint) // to allow access from python
+			(initContactOnPositiveFacetSide)
 			// hasShear
-			REGISTER_ATTRIBUTE(hasShear);
-			REGISTER_ATTRIBUTE(pos1);
-			REGISTER_ATTRIBUTE(pos2);
-			REGISTER_ATTRIBUTE(ori1);
-			REGISTER_ATTRIBUTE(ori2);
-			REGISTER_ATTRIBUTE(cp1rel);
-			REGISTER_ATTRIBUTE(cp2rel);
-			REGISTER_ATTRIBUTE(d1);
-			REGISTER_ATTRIBUTE(d2);
-			REGISTER_ATTRIBUTE(d0);
-			REGISTER_ATTRIBUTE(initRelOri12);
-		}
-	REGISTER_CLASS_NAME(SpheresContactGeometry);
-	REGISTER_BASE_CLASS_NAME(InteractionGeometry);
+			(hasShear)
+			(pos1)
+			(pos2)
+			(ori1)
+			(ori2)
+			(cp1rel)
+			(cp2rel)
+			(d1)
+			(d2)
+			(d0)
+			(initRelOri12));
+	REGISTER_CLASS_AND_BASE(SpheresContactGeometry,InteractionGeometry);
 	REGISTER_CLASS_INDEX(SpheresContactGeometry,InteractionGeometry);
 };
 

Modified: trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/pkg/dem/Engine/EngineUnit/InteractingFacet2InteractingSphere4SpheresContactGeometry.cpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -41,7 +41,7 @@
 	Matrix3r facetAxisT; se31.orientation.ToRotationMatrix(facetAxisT); 
 	Matrix3r facetAxis = facetAxisT.Transpose();
 	// local orientation
-	Vector3r cl = facetAxis*(se32.position - se31.position); 
+	Vector3r cl = facetAxis*(se32.position - se31.position);  // "contact line" in facet-local coords
 
 	//
 	// BEGIN everything in facet-local coordinates
@@ -49,10 +49,17 @@
 
 	Vector3r normal = facet->nf;
 	Real L = normal.Dot(cl);
-	if (L<0) { normal=-normal; L=-L; }
-	Real sphereRadius = static_cast<InteractingSphere*>(cm2.get())->radius;
 
+	bool newOnPositiveSide; // temp to save what will be maybe needed for new contact
+	if(c->interactionGeometry){ // contact already exists, use old data here
+		const bool& icopfs=YADE_CAST<SpheresContactGeometry*>(c->interactionGeometry.get())->initContactOnPositiveFacetSide;
+		if(!icopfs) { normal=-normal; L=-L; } // original contact on the other side, reverse the sense now
+	} else {
+		if (L<0) { normal=-normal; L=-L; newOnPositiveSide=false;} // new contact on the negative face, reverse and save that information so that since now this contact is always reversed
+		else newOnPositiveSide=true;
+	}
 
+	Real sphereRadius = static_cast<InteractingSphere*>(cm2.get())->radius;
 	if (L > sphereRadius)  return false; // no contact
 
 	Vector3r cp = cl - L*normal;
@@ -108,7 +115,8 @@
 		if (c->interactionGeometry)
 			scm = YADE_PTR_CAST<SpheresContactGeometry>(c->interactionGeometry);
 		else
-			scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
+			{ scm = shared_ptr<SpheresContactGeometry>(new SpheresContactGeometry());
+			scm->initContactOnPositiveFacetSide=newOnPositiveSide; }
 	  
 		normal = facetAxisT*normal; // in global orientation
 		scm->contactPoint = se32.position - (sphereRadius-0.5*penetrationDepth)*normal; 

Modified: trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp
===================================================================
--- trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/pkg/dem/RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp	2009-01-20 18:30:38 UTC (rev 1630)
@@ -79,12 +79,17 @@
 
 	if(sc->hasShear){
 		Vector3r pos1=sc->pos1, pos2=sc->pos2, contPt=sc->contPt();
-		//Vector3r contPt=se31.position+(sc->d1/sc->d0)*(se32.position-se31.position); // must be recalculated to not be unscaled if scaling displacements ...
-		GLUtils::GLDrawLine(pos1,pos2,Vector3r(.5,.5,.5));
-		Vector3r bend; Real tors;
-		sc->bendingTorsionRel(bend,tors);
-		GLUtils::GLDrawLine(contPt,contPt+10*sc->radius1*(bend+sc->normal*tors),Vector3r(1,0,0));
 		#if 0
+			GLUtils::GLDrawArrow(contPt,contPt+sc->normal*.5*sc->d0); // normal of the contact
+		#endif
+		#if 0
+			//Vector3r contPt=se31.position+(sc->d1/sc->d0)*(se32.position-se31.position); // must be recalculated to not be unscaled if scaling displacements ...
+			GLUtils::GLDrawLine(pos1,pos2,Vector3r(.5,.5,.5));
+			Vector3r bend; Real tors;
+			sc->bendingTorsionRel(bend,tors);
+			GLUtils::GLDrawLine(contPt,contPt+10*sc->radius1*(bend+sc->normal*tors),Vector3r(1,0,0));
+		#endif
+		#if 0
 		GLUtils::GLDrawNum(bend[0],contPt-.2*sc->normal*sc->radius1,Vector3r(1,0,0));
 		GLUtils::GLDrawNum(bend[1],contPt,Vector3r(0,1,0));
 		GLUtils::GLDrawNum(bend[2],contPt+.2*sc->normal*sc->radius1,Vector3r(0,0,1));
@@ -96,14 +101,16 @@
 		//TRVAR4(pos1,sc->ori1,pos2,sc->ori2);
 		//TRVAR2(sc->cp2rel,pos2+(sc->ori2*sc->cp2rel*Vector3r::UNIT_X*sc->d2));
 		// contact point to projected points
-		Vector3r ptTg1=sc->contPtInTgPlane1(), ptTg2=sc->contPtInTgPlane2();
-		TRVAR3(ptTg1,ptTg2,sc->normal)
-		//GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1)); GLUtils::GLDrawLine(pos1,contPt+ptTg1,Vector3r(0,.5,1));
-		//GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5)); GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
-		// projected shear
-		GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
-		// 
-		//GLUtils::GLDrawNum(sc->epsN(),contPt,Vector3r(1,1,1));
+		#if 0
+			Vector3r ptTg1=sc->contPtInTgPlane1(), ptTg2=sc->contPtInTgPlane2();
+			TRVAR3(ptTg1,ptTg2,sc->normal)
+			//GLUtils::GLDrawLine(contPt,contPt+ptTg1,Vector3r(0,.5,1)); GLUtils::GLDrawLine(pos1,contPt+ptTg1,Vector3r(0,.5,1));
+			//GLUtils::GLDrawLine(contPt,contPt+ptTg2,Vector3r(0,1,.5)); GLUtils::GLDrawLine(pos2,contPt+ptTg2,Vector3r(0,1,.5));
+			// projected shear
+			GLUtils::GLDrawLine(contPt+ptTg1,contPt+ptTg2,Vector3r(1,1,1));
+			// 
+			//GLUtils::GLDrawNum(sc->epsN(),contPt,Vector3r(1,1,1));
+		#endif
 	}
 
 

Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/pkg/dem/SConscript	2009-01-20 18:30:38 UTC (rev 1630)
@@ -385,7 +385,7 @@
 		['RenderingEngine/GLDrawSpheresContactGeometry/GLDrawSpheresContactGeometry.cpp'],
 		LIBS=env['LIBS']+['SpheresContactGeometry',
 			'ElasticContactInteraction',
-			'yade-opengl']),
+			'yade-opengl','$QGLVIEWER_LIB']),
 
 	env.SharedLibrary('GLDrawSDECLinkGeometry',
 		['RenderingEngine/GLDrawSDECLinkGeometry/GLDrawSDECLinkGeometry.cpp'],

Added: trunk/scripts/test-sphere-facet.py
===================================================================
--- trunk/scripts/test-sphere-facet.py	2009-01-19 15:19:01 UTC (rev 1629)
+++ trunk/scripts/test-sphere-facet.py	2009-01-20 18:30:38 UTC (rev 1630)
@@ -0,0 +1,65 @@
+# © Václav Šmilauer <eudoxos@xxxxxxxx>
+#
+# Test case for sphere-facet interaction preserving the original contact orientation.
+# Z-gravity is being increased every 4000 iterations, the sphere dives more into the facet and stabilizes,
+# etc. This process continues even if the sphere center passes on the other side of the facet and,
+# (if distant transient interactions are allowed in the collider) if the sphere passes in its entirety to
+# the other side of the facet. The interaction, however, still pushes in the same sense.
+#
+# After the gravity reaches some value, it is reset and the sphere should be pushed from the facet towards
+# its original position. When the contact on the original side is lost, the interaction should be deleted.
+#
+#
+# The only tunable sign places the sphere either on the top ot at the bottom of the facet
+# and sets gravity accordingly. It can be +1 or -1
+#
+sign=-1
+#
+
+O.initializers=[
+	StandAloneEngine('PhysicalActionContainerInitializer'),
+	MetaEngine('BoundingVolumeMetaEngine',[EngineUnit('InteractingSphere2AABB'),EngineUnit('InteractingBox2AABB'),EngineUnit('MetaInteractingGeometry2AABB')])
+	]
+O.engines=[
+	StandAloneEngine('PhysicalActionContainerReseter'),
+	MetaEngine('BoundingVolumeMetaEngine',[
+		EngineUnit('InteractingSphere2AABB'),
+		EngineUnit('InteractingBox2AABB'),
+		EngineUnit('MetaInteractingGeometry2AABB')
+	]),
+	StandAloneEngine('PersistentSAPCollider',{'haveDistantTransient':True}),
+	MetaEngine('InteractionGeometryMetaEngine',[
+		EngineUnit('InteractingSphere2InteractingSphere4SpheresContactGeometry',{'hasShear':True,'interactionDetectionFactor':1.4}),
+		EngineUnit('InteractingFacet2InteractingSphere4SpheresContactGeometry',{'hasShear':True}),
+	]),
+	MetaEngine('InteractionPhysicsMetaEngine',[EngineUnit('SimpleElasticRelationships')]),
+	StandAloneEngine('ElasticContactLaw'),
+	DeusExMachina('GravityEngine',{'gravity':[0,0,-sign*500],'label':'gravitator'}),
+	DeusExMachina("NewtonsDampedLaw",{'damping':0.8}),
+	StandAloneEngine('PeriodicPythonRunner',{'iterPeriod':4000,'command':'setGravity()'}),
+	]
+O.bodies.append([
+	utils.facet([[-1,-1,0],[1,-1,0],[0,1,0]],dynamic=False,color=[1,0,0],young=1e3),
+	utils.sphere([0,0,sign*.49999],radius=.5,young=1e3,wire=True,density=1),
+])
+O.miscParams=[Generic('GLDrawSphere',{'glutUse':True})]
+O.saveTmp('init')
+O.dt=1e-4
+
+def setGravity():
+	gz=gravitator["gravity"][2]
+	gravitator["gravity"]=[0,0,1.05*gz]
+	if abs(gz)>=2500:
+		print "Gravity reset & slow down"
+		O.dt=1e-6;
+		gravitator["gravity"]=[0,0,0]
+	if abs(gz)>0: print gz
+
+try:
+	from yade import qt
+	renderer=qt.Renderer()
+	renderer['Interaction_geometry']=True
+	qt.Controller()
+except ImportError: pass
+
+