yade-dev team mailing list archive

[noreply@xxxxxxxxxxxxx]

------------------------------------------------------------
revno: 2474
committer: Václav Šmilauer <eu@xxxxxxxx>
branch nick: yade
timestamp: Tue 2010-10-12 14:28:40 +0200
message:
  1. Fix QT4CXX key error, thanks to Anton for reporting.
modified:
  gui/SConscript
  pkg/dem/Engine/GlobalEngine/DomainLimiter.cpp
  pkg/dem/Engine/GlobalEngine/DomainLimiter.hpp
  scripts/test/law-test.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 'gui/SConscript'
--- gui/SConscript	2010-10-12 11:47:16 +0000
+++ gui/SConscript	2010-10-12 12:28:40 +0000
@@ -2,6 +2,8 @@
 Import('*')
 linkPlugins=env['linkPlugins']
 
+haveQt4Cxx=('QT4CXX' in env and env['QT4CXX'])
+
 if 'qt4' in env['features']:
 	env.Install('$LIBDIR/py/yade/qt',[
 		env.File('qt4/img_rc.py'),
@@ -9,7 +11,7 @@
 		env.File('qt4/SerializableEditor.py'),
 		env.File('qt4/Inspector.py'),
 		env.File('qt4/__init__.py'),
-		env.SharedLibrary('_GLViewer',['qt4/GLViewer.cpp','qt4/_GLViewer.cpp','qt4/OpenGLManager.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+[env['QGLVIEWER_LIB']]+linkPlugins(['PeriodicEngines']),RPATH=env['RPATH']+[env.Literal('\\$$ORIGIN/../../../gui')],CXX=env['QT4CXX'] if env['QT4CXX'] else env['CXX'],CXXFLAGS=[f for f in env['CXXFLAGS'] if not f.startswith('-Q')])
+		env.SharedLibrary('_GLViewer',['qt4/GLViewer.cpp','qt4/_GLViewer.cpp','qt4/OpenGLManager.cpp'],SHLIBPREFIX='',LIBS=env['LIBS']+[env['QGLVIEWER_LIB']]+linkPlugins(['PeriodicEngines']),RPATH=env['RPATH']+[env.Literal('\\$$ORIGIN/../../../gui')],CXX=env['QT4CXX'] if haveQt4Cxx else env['CXX'],CXXFLAGS=[f for f in env['CXXFLAGS'] if not f.startswith('-Q')] if haveQt4Cxx else env['CXXFLAGS'])
 	])
 	env.Command('qt4/img_rc.py','qt4/img.qrc','pyrcc4 -o $buildDir/gui/qt4/img_rc.py gui/qt4/img.qrc')
 	env.Command('qt4/ui_controller.py','qt4/controller.ui','pyuic4 -o $buildDir/gui/qt4/ui_controller.py gui/qt4/controller.ui')

=== modified file 'pkg/dem/Engine/GlobalEngine/DomainLimiter.cpp'
--- pkg/dem/Engine/GlobalEngine/DomainLimiter.cpp	2010-10-12 08:33:10 +0000
+++ pkg/dem/Engine/GlobalEngine/DomainLimiter.cpp	2010-10-12 12:28:40 +0000
@@ -89,7 +89,7 @@
 		if(scGeom){
 			scGeom->rotate(axY); scGeom->rotate(axZ);
 			scGeom->rotate(shearTot);
-			shearTot-=scGeom->shearIncrement();
+			shearTot+=scGeom->shearIncrement();
 			ptGeom=Vector3r(-scGeom->penetrationDepth,shearTot.dot(axY),shearTot.dot(axZ));
 		}
 		else{ // d3dGeom
@@ -103,30 +103,57 @@
 	trsf.row(0)=axX; trsf.row(1)=axY; trsf.row(2)=axZ;
 	trsfQ=Quaternionr(trsf);
 	contPt=gsc->contactPoint;
+	Real refLength=gsc->refR1+gsc->refR2;
 	
 	// here we go ahead, finally
 	Vector3r val=linearInterpolate<Vector3r,int>(step,_pathT,_pathV,_interpPos);
 	Vector3r valDiff=val-prevVal; prevVal=val;
-	if(pathIsRel){
+	if(displIsRel){
 		LOG_DEBUG("Relative diff is "<<valDiff<<" (will be normalized by "<<gsc->refR1+gsc->refR2<<")");
-		valDiff*=gsc->refR1+gsc->refR2;
+		valDiff*=refLength;
 	}
 	LOG_DEBUG("Absolute diff is "<<valDiff);
 	ptOurs+=valDiff;
 
-	#if 0
-		// we compute displacement for the 2nd particles; cff1 and cff2 determine how to apply it on both particles, depending on sense:
-		// 0 is the symmetric case (-.5,.5); -1 applies only to the first particle, i.e. (-1,0), +1 only to the second one (0,+1)
-		Real r1_r2=gsc->radius1/gsc->radius2; // ratio of radii
-		Real cff1=(sense<0?-1:(sense>0?0:-.5)), cff2=(sense<0?0:(sense<0?0:.5));
-	#endif
 	// shear is applied as rotation of id2: dε=r₁dθ → dθ=dε/r₁;
-	Real rot2y=valDiff[1]/gsc->refR2, rot2z=valDiff[2]/gsc->refR2;
-	Vector3r angVel2=(rot2y*axY+rot2z*axZ)/scene->dt;
-	Vector3r linVel2=axX*valDiff[0]/scene->dt;
-	state2->angVel=angVel2;
-	state2->vel=linVel2;
-	LOG_DEBUG("Body #"<<id2<<", setting vel="<<linVel2<<", angVel="<<angVel2);
+	Vector3r vel[2],angVel[2];
+	for(int i=0; i<2; i++){
+		int sign=(i==0?-1:1);
+		Real weight=(i==0?1-idWeight:idWeight);
+		Real radius=(i==0?gsc->refR1:gsc->refR2);
+		Vector3r diff=sign*valDiff*weight;
+
+		// normal displacement
+		vel[i]=axX*diff[0]/scene->dt;
+		// shear rotation: angle
+		Real rotZ=rotWeight*diff[1]/radius, rotY=rotWeight*diff[2]/radius;
+		angVel[i]=(rotY*axY+rotZ*axZ)/scene->dt;
+		// shear displacement
+		//Real arcAngleY=atan(weight*(1-rotWeight)*sign*valDiff[1]/radius), arcAngleZ=atan(weight*(1-rotWeight)*sign*valDiff[2]/radius);
+		Real arcAngleY=(1-rotWeight)*diff[1]/radius, arcAngleZ=(1-rotWeight)*diff[2]/radius;
+		//Real arcAngleY=arcLenY/radius, arcAngleZ=arcLenZ/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
+		vel[i]-=sign*axX*.5*radius*((1-cos(arcAngleY))+(1-cos(arcAngleZ)))/scene->dt;
+
+		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];
+	LOG_DEBUG("Body #"<<id1<<", setting vel="<<vel[0]<<", angVel="<<angVel[0]);
+	LOG_DEBUG("Body #"<<id2<<", setting vel="<<vel[1]<<", angVel="<<angVel[1]);
+
+	#if 0
+		Real rot2y=valDiff[1]/gsc->refR2, rot2z=valDiff[2]/gsc->refR2;
+		Vector3r angVel2=(rot2y*axY+rot2z*axZ)/scene->dt;
+		Vector3r linVel2=axX*valDiff[0]/scene->dt;
+		state2->angVel=angVel2;
+		state2->vel=linVel2;
+		LOG_DEBUG("Body #"<<id2<<", setting vel="<<linVel2<<", angVel="<<angVel2);
+	#endif
 
 	step++;
 }
@@ -143,26 +170,34 @@
 	if(!tester){ FOREACH(shared_ptr<Engine> e, scene->engines){ tester=dynamic_pointer_cast<LawTester>(e); if(tester) break; } }
 	if(!tester){ LOG_ERROR("No LawTester in O.engines, killing myself."); dead=true; return; }
 	//if(tester->renderLength<=0) return;
-	glColor3v(Vector3r(0,1,0));
+	glColor3v(Vector3r(1,0,1));
+
+	glLineWidth(2.0f);
+	glEnable(GL_LINE_SMOOTH);
+	glLoadIdentity();
+
+	glTranslatev(tester->contPt);
+
 	glBegin(GL_LINES); glVertex3v(Vector3r::Zero()); glVertex3v(.1*Vector3r::Ones()); glEnd(); 
-	//cerr<<".";
-	return; 
+	GLUtils::GLDrawText(string("This is the contact point!"),Vector3r::Zero(),Vector3r(1,0,1));
+
+	cerr<<".";
 
 	glBegin(GL_LINES);
-		glVertex3v(tester->contPt);
-		glVertex3v(tester->contPt+tester->renderLength*tester->axX);
-		glVertex3v(tester->contPt);
-		glVertex3v(tester->contPt+tester->renderLength*tester->axY);
-		glVertex3v(tester->contPt);
-		glVertex3v(tester->contPt+tester->renderLength*tester->axZ);
-		glVertex3v(tester->contPt);
-		glVertex3v(tester->contPt+tester->ptOurs);
+		glVertex3v(Vector3r::Zero());
+		glVertex3v(tester->renderLength*tester->axX);
+		glVertex3v(Vector3r::Zero());
+		glVertex3v(tester->renderLength*tester->axY);
+		glVertex3v(Vector3r::Zero());
+		glVertex3v(tester->renderLength*tester->axZ);
+		glVertex3v(Vector3r::Zero());
+		glVertex3v(tester->ptOurs);
 	glEnd();
-	#if 0
-	GLUtils::GLDrawArrow(tester->contPt,tester->contPt+tester->renderLength*tester->axX,Vector3r(1,0,0));
-	GLUtils::GLDrawArrow(tester->contPt,tester->contPt+tester->renderLength*tester->axY,Vector3r(0,1,0));
-	GLUtils::GLDrawArrow(tester->contPt,tester->contPt+tester->renderLength*tester->axZ,Vector3r(0,0,1));
-	GLUtils::GLDrawArrow(tester->contPt,tester->contPt+tester->ptOurs,Vector3r(1,1,0));
+	#if 1
+	GLUtils::GLDrawArrow(Vector3r::Zero(),tester->contPt+tester->renderLength*tester->axX,Vector3r(1,0,0));
+	GLUtils::GLDrawArrow(Vector3r::Zero(),tester->contPt+tester->renderLength*tester->axY,Vector3r(0,1,0));
+	GLUtils::GLDrawArrow(Vector3r::Zero(),tester->contPt+tester->renderLength*tester->axZ,Vector3r(0,0,1));
+	GLUtils::GLDrawArrow(Vector3r::Zero(),tester->contPt+tester->ptOurs,Vector3r(1,1,0));
 	GLUtils::GLDrawArrow(Vector3r::Zero(),Vector3r::Ones(),Vector3r(.5,.5,1));
 	#endif
 }

=== modified file 'pkg/dem/Engine/GlobalEngine/DomainLimiter.hpp'
--- pkg/dem/Engine/GlobalEngine/DomainLimiter.hpp	2010-10-01 15:38:55 +0000
+++ pkg/dem/Engine/GlobalEngine/DomainLimiter.hpp	2010-10-12 12:28:40 +0000
@@ -25,7 +25,7 @@
 		((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|"))
-		((bool,pathIsRel,true,,"Whether values in *path* are normalized by contact length."))
+		((bool,displIsRel,true,,"Whether displacement values in *path* are normalized by reference contact length (r1+r2 for 2 spheres)."))
 		((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."))
@@ -41,9 +41,8 @@
 		((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>`."))
 		((Real,renderLength,0,,"Characteristic length for the purposes of rendering, set equal to the smaller radius."))
 		((Vector3r,contPt,Vector3r::Zero(),,"Contact point (for rendering only)"))
-		//((int,sense,1,,"Determines what particle is moved: negative for id1, positive for id2, 0 for both."))
-		// applyWeight
-		// shearRotWeight
+		((Real,idWeight,1,,"Float ∈〈0,1〉 determining on which particle are displacements applied (0 for id1, 1 for id2); intermediate values will apply respective part to each of them."))
+		((Real,rotWeight,1,,"Float ∈〈0,1〉 determining whether shear displacement is applied as rotation or displacement on arc (0 is displacemetn-only, 1 is rotation-only)."))
 		// reset force components along individual axes, instead of blocking DOFs which have no specific direction (for the force control)
 	);
 };

=== modified file 'scripts/test/law-test.py'
--- scripts/test/law-test.py	2010-09-30 18:00:41 +0000
+++ scripts/test/law-test.py	2010-10-12 12:28:40 +0000
@@ -13,14 +13,15 @@
 #
 from yade import utils,plot
 import random, yade.log
-#yade.log.setLevel('LawTester',yade.log.TRACE)
+yade.log.setLevel('LawTester',yade.log.TRACE)
 
 # sphere's radii
 r1,r2=.1,.2
 # place sphere 1 at the origin
 pt1=Vector3(0,0,0)
 # random orientation of the interaction
-normal=Vector3(random.random()-.5,random.random()-.5,random.random()-.5)
+# normal=Vector3(random.random()-.5,random.random()-.5,random.random()-.5)
+normal=Vector3(1,0,0)
 O.bodies.append([
 	utils.sphere(pt1,r1,wire=True),
 	utils.sphere(pt1+.999999*(r1+r2)*normal.normalized(),r2,wire=True)
@@ -35,9 +36,10 @@
 		[Law2_ScGeom_FrictPhys_Basic()]
 	),
 	LawTester(ids=[0,1],path=[
+		#(-.1,0,0),(-.1,.1,0) #,(-.1,0,.1)
 		(-.1,0,0),(-.1,.1,0),(0,.1,0), # towards, shear, back to intial normal distance
 		(-.02,.1,.1),(-.02,-.1,.1),(-.02,-.1,-.1),(-.02,.1,-.1),(-.02,.1,.1) # go in square in the shear plane without changing normal deformation
-		],pathSteps=[100],doneHook='tester.dead=True; O.pause()',label='tester'),
+		],pathSteps=[10],doneHook='tester.dead=True; O.pause()',label='tester',rotWeight=1,idWeight=1),
 	PyRunner(iterPeriod=1,command='addPlotData()'),
 	NewtonIntegrator()
 ]
@@ -45,7 +47,7 @@
 def addPlotData():
 	i=O.interactions[0,1]
 	plot.addData(un=tester.ptOurs[0],us1=tester.ptOurs[1],us2=tester.ptOurs[2],ung=tester.ptGeom[0],us1g=tester.ptGeom[1],us2g=tester.ptGeom[2],i=O.iter,Fs=i.phys.shearForce.norm(),Fn=i.phys.normalForce.norm())
-plot.plots={'us1':('us2',),'Fn':('Fs',),'i':('un','us1','us2'),' i':('Fs','Fn')}  #'ung','us1g','us2g'
+plot.plots={'us1':('us2',),'Fn':('Fs',),'i':('un','us1','us2'),' i':('Fs','Fn'),'  i':('ung','us1g','us2g')}  #'ung','us1g','us2g'
 plot.plot(subPlots=True)
 
 try:
@@ -54,6 +56,7 @@
 	rr=qt.Renderer()
 	rr.extraDrawers=[GlExtra_LawTester()]
 except ImportError: pass
+O.dt=1e2
 
 O.saveTmp()
 O.run()