yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2008-12-01 20:47:18 +0100 (Mon, 01 Dec 2008)
New Revision: 1587

Added:
   trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp
   trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp
   trunk/scripts/test/
   trunk/scripts/test/Se3Interpolator.py
Modified:
   trunk/extra/Brefcom.cpp
   trunk/extra/usct/UniaxialStrainControlledTest.cpp
   trunk/extra/usct/UniaxialStrainControlledTest.hpp
   trunk/gui/py/PythonUI_rc.py
   trunk/pkg/common/Engine/StandAloneEngine/PeriodicEngines.hpp
   trunk/pkg/common/Engine/StandAloneEngine/ResetPositionEngine.cpp
   trunk/pkg/common/SConscript
Log:
1. Added Se3Interpolator engine, with scripts/test/Se3Interpolator.py
2. Added from math import * to PythonUI_rc.py so that math functions are readily accessible everywhere
3. Implemented stopStrain in UniaxialStrainControlledTest
4. other small fixes


Modified: trunk/extra/Brefcom.cpp
===================================================================
--- trunk/extra/Brefcom.cpp	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/extra/Brefcom.cpp	2008-12-01 19:47:18 UTC (rev 1587)
@@ -156,6 +156,7 @@
 
 		epsN=contGeom->epsN();
 		epsT=contGeom->epsT();
+		//if(I->getId1()==0 && (Omega::instance().getCurrentIteration()%10==0)) LOG_INFO("##"<<I->getId1()<<"+"<<I->getId2()<<": "<<" dist2-1="<<(contGeom->pos1-contGeom->pos2).Length()-1.<<", d0="<<contGeom->d0<<", epsN="<<epsN<<", |epsT|="<<epsT.Length());
 
 		#ifdef BREFCOM_MATERIAL_MODEL
 			BREFCOM_MATERIAL_MODEL

Modified: trunk/extra/usct/UniaxialStrainControlledTest.cpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.cpp	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/extra/usct/UniaxialStrainControlledTest.cpp	2008-12-01 19:47:18 UTC (rev 1587)
@@ -84,8 +84,19 @@
 	if(posIds.size()==0 || negIds.size()==0) return;
 	// linearly increase strain to the desired value
 	if(abs(currentStrainRate)<abs(strainRate))currentStrainRate+=strainRate*.005; else currentStrainRate=strainRate;
-	// how much do we move; in the symmetric case, half of the strain is applied at each end;
-	Real dAX=(asymmetry==0?.5:1)*currentStrainRate*originalLength*Omega::instance().getTimeStep();
+	// how much do we move (in total, symmetry handled below)
+	Real dAX=currentStrainRate*originalLength*Omega::instance().getTimeStep();
+	if(!isnan(stopStrain)){
+		Real axialLength=axisCoord(posIds[0])-axisCoord(negIds[0]);
+		Real newStrain=(axialLength+dAX)/originalLength-1;
+		if((newStrain*stopStrain>0) && abs(newStrain)>=stopStrain){ // same sign of newStrain and stopStrain && over the limit from below in abs values
+			dAX=originalLength*(stopStrain+1)-axialLength;
+			LOG_INFO("Reached stopStrain "<<stopStrain<<", deactivating self and stopping in "<<idleIterations+1<<" iterations.");
+			this->active=false;
+			rootBody->stopAtIteration=Omega::instance().getCurrentIteration()+1+idleIterations;
+		}
+	}
+	if(asymmetry==0) dAX*=.5; // apply half on both sides if straining symetrically
 	for(size_t i=0; i<negIds.size(); i++){
 		if(asymmetry==0 || asymmetry==-1 /* for +1, don't move*/) negCoords[i]-=dAX;
 		axisCoord(negIds[i])=negCoords[i]; // update current position
@@ -102,9 +113,6 @@
 	// reverse if we're over the limit strain
 	if(notYetReversed && limitStrain!=0 && ((currentStrainRate>0 && strain>limitStrain) || (currentStrainRate<0 && strain<limitStrain))) { currentStrainRate*=-1; notYetReversed=false; LOG_INFO("Reversed strain rate to "<<currentStrainRate); }
 
-	// if(strain==stopAtStrain){LOG_INFO("Prescibed strain "<<stopAtStrain<<" reached, pausing."); }
-	if(!isnan(stopAtStrain) && ((strain*stopAtStrain>0) && abs(strain)>abs(stopAtStrain))) { strain=stopAtStrain;	rootBody->stopAtIteration=Omega::instance().getCurrentIteration()+2; LOG_INFO("Prescribed strain "<<stopAtStrain<<" reached, will stop in 2 iterations"); }
-
 	if(Omega::instance().getCurrentIteration()%10==0 ) {
 		computeAxialForce(rootBody);
 		#if 0
@@ -297,12 +305,6 @@
 
 	rootBody->engines.push_back(shared_ptr<BrefcomDamageColorizer>(new BrefcomDamageColorizer));
 
-	shared_ptr<PositionOrientationRecorder> por(new PositionOrientationRecorder);
-	por->outputFile="/tmp/usct-traction";
-	por->interval=300;
-	por->saveRgb=true;
-	rootBody->engines.push_back(por);
-
 }
 
 

Modified: trunk/extra/usct/UniaxialStrainControlledTest.hpp
===================================================================
--- trunk/extra/usct/UniaxialStrainControlledTest.hpp	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/extra/usct/UniaxialStrainControlledTest.hpp	2008-12-01 19:47:18 UTC (rev 1587)
@@ -80,9 +80,10 @@
 		Real& axisCoord(body_id_t id){ return Body::byId(id)->physicalParameters->se3.position[axis]; };
 		void init();
 	public:
+		virtual bool isActivated(){return active;}
 		Real strainRate,currentStrainRate;
 		//! strain at which we will pause simulation; inactive (nan) by default; must be reached from below (in absolute value)
-		Real stopAtStrain;
+		Real stopStrain;
 		//! distance of reference bodies in the direction of axis before straining started
 		Real originalLength;
 		vector<Real> originalWidths;
@@ -100,6 +101,10 @@
 		bool blockDisplacements;
 		//! Whether rotations of boundary bodies are blocked.
 		bool blockRotations;
+		//! Are we activated?
+		bool active;
+		//! Number of iterations that will pass without straining activity after stopStrain has been reached (default: 0)
+		long idleIterations;
 
 		/** bodies on which straining will be applied (on the positive and negative side of axis) */
 		vector<body_id_t> posIds, negIds;
@@ -126,11 +131,13 @@
 		Real strain, avgStress;
 
 		virtual void applyCondition(MetaBody* rootBody);
-		UniaxialStrainer(){axis=2; asymmetry=0; currentStrainRate=0; originalLength=-1; limitStrain=0; notYetReversed=true; crossSectionArea=-1; needsInit=true; /* sensorsPusher=shared_ptr<UniaxialStrainSensorPusher>(); */ recordFile="/tmp/usct.data"; strain=avgStress=/*avgTransStrain=*/0; blockRotations=false; blockDisplacements=false;  stopAtStrain=numeric_limits<Real>::quiet_NaN();};
+		UniaxialStrainer(){axis=2; asymmetry=0; currentStrainRate=0; originalLength=-1; limitStrain=0; notYetReversed=true; crossSectionArea=-1; needsInit=true; /* sensorsPusher=shared_ptr<UniaxialStrainSensorPusher>(); */ recordFile="/tmp/usct.data"; strain=avgStress=/*avgTransStrain=*/0; blockRotations=false; blockDisplacements=false;  stopStrain=numeric_limits<Real>::quiet_NaN(); active=true; idleIterations=0; };
 		virtual ~UniaxialStrainer(){};
 		REGISTER_ATTRIBUTES_WITH_BASE(DeusExMachina,
 				(strainRate) 
-				(stopAtStrain) 
+				(stopStrain) 
+				(active)
+				(idleIterations)
 				(currentStrainRate) 
 				(axis) 
 				(asymmetry) 

Modified: trunk/gui/py/PythonUI_rc.py
===================================================================
--- trunk/gui/py/PythonUI_rc.py	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/gui/py/PythonUI_rc.py	2008-12-01 19:47:18 UTC (rev 1587)
@@ -7,6 +7,7 @@
 sys.excepthook=sys.__excepthook__ # apport on ubuntu overrides this, we don't need it
 # sys.path.insert(0,runtime.prefix+'/lib/yade'+runtime.suffix+'/extra')
 
+from math import *
 from yade.wrapper import *
 from yade import runtime
 from yade import utils

Added: trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.cpp	2008-12-01 19:47:18 UTC (rev 1587)
@@ -0,0 +1,37 @@
+
+#include"Se3Interpolator.hpp"
+#include<yade/pkg-common/PeriodicEngines.hpp>
+
+YADE_PLUGIN("Se3Interpolator")
+CREATE_LOGGER(Se3Interpolator);
+
+void Se3Interpolator::applyCondition(MetaBody* mb){
+	assert(subscribedBodies.size()==1);
+	const shared_ptr<Body>&b=Body::byId(subscribedBodies[0],mb);
+	Omega& O=Omega::instance();
+	if(!started){
+		started=true;
+		start=b->physicalParameters->se3;
+		if(rotRelative) goal.orientation=goal.orientation*start.orientation;
+		int numGoalEndpoints=0;
+		if(goalVirt>0){startVirt=O.getSimulationTime(); numGoalEndpoints++;}
+		if(goalReal>0){startReal=PeriodicEngine::getClock(); numGoalEndpoints++;}
+		if(goalIter>0){startIter=O.getCurrentIteration(); numGoalEndpoints++;}
+		assert(numGoalEndpoints==1);
+		LOG_DEBUG("Init'ed.");
+	}
+	Real t=-1; // t∈[0,1]
+	if(goalVirt>0) t=(O.getSimulationTime()-startVirt)/(goalVirt-startVirt);
+	if(goalReal>0) t=(PeriodicEngine::getClock()-startReal)/(goalReal-startReal);
+	if(goalIter>0) t=(O.getCurrentIteration()-startIter)/(goalIter-startIter);
+	assert(t>=0);
+	LOG_DEBUG("t="<<t);
+	t=min(t,1.);
+	b->physicalParameters->se3.position=start.position+t*(goal.position-start.position);
+	b->physicalParameters->se3.orientation=Quaternionr().Slerp(t,start.orientation,goal.orientation);
+	if(t>=1.){
+		done=true;
+		LOG_DEBUG("Goal reached.");
+		if(!goalHook.empty()){ PyGILState_STATE gstate; gstate=PyGILState_Ensure(); PyRun_SimpleString(goalHook.c_str()); PyGILState_Release(gstate); }
+	}
+}

Added: trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp
===================================================================
--- trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/pkg/common/Engine/DeusExMachina/Se3Interpolator.hpp	2008-12-01 19:47:18 UTC (rev 1587)
@@ -0,0 +1,60 @@
+// 2008 © Václav Šmilauer <eudoxos@xxxxxxxx> 
+#pragma once
+
+#include<yade/core/DeusExMachina.hpp>
+#include<yade/core/PhysicalParameters.hpp>
+#include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<Python.h>
+
+/* Engine interpolating between starting (current) and goal (given) se3, both position and orientation.
+ *
+ * The interpolation is linear, either in iteration space or (real/virtual) time space.
+ *
+ * A given python code can be executed when the ending se3 is reached and the engine deactivates itself.
+ *
+ * Note that eactly one body must be subscribed to this engine, since interpolating positions
+ * on multiple bodies is meaningless.
+ *
+ * Attributes that should be assigned:
+ * 	subscribedBodies (exactly one)
+ * 	goal (7-tuple of position and orientation)
+ * 	startVirt XOR startReal XOR startIter (at what time point to start)
+ * 	goalVirt XOR goalReal XOR goalIter (at what time point to reach desired se3)
+ * 	goalHook (optional: python command to be run when finished)
+ * 	rotRelative (optional: rotation part of goal is relative to the rotation at the beginning)
+ */
+class Se3Interpolator: public DeusExMachina {
+	public:
+		bool done,started;
+		Se3r start;
+		Se3r goal;
+		bool rotRelative;
+		Real goalVirt, goalReal; long goalIter;
+		long startVirt, startReal; long startIter;
+		string goalHook;
+		virtual bool isActivated(){return !done;}
+		Se3Interpolator(): done(false), started(false), rotRelative(false), goalVirt(0), goalReal(0), goalIter(0) {};
+		virtual ~Se3Interpolator(){};
+		virtual void applyCondition(MetaBody* mb);
+	protected:
+		REGISTER_ATTRIBUTES_WITH_BASE(DeusExMachina,
+			(done)
+			(started)
+			(start)
+			(goal)
+			(rotRelative)
+			(goalVirt)
+			(goalReal)
+			(goalIter)
+			(startVirt)
+			(startReal)
+			(startIter)
+			(goalHook)
+		);
+		REGISTER_CLASS_AND_BASE(Se3Interpolator,DeusExMachina);
+		DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Se3Interpolator);
+	
+
+

Modified: trunk/pkg/common/Engine/StandAloneEngine/PeriodicEngines.hpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/PeriodicEngines.hpp	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/pkg/common/Engine/StandAloneEngine/PeriodicEngines.hpp	2008-12-01 19:47:18 UTC (rev 1587)
@@ -11,9 +11,8 @@
  * the number of activations was already reached, no action will be called even if any of active period has elapsed.
  */
 class PeriodicEngine:  public StandAloneEngine {
-	protected:
-		static Real getClock(){ timeval tp; gettimeofday(&tp,NULL); return tp.tv_sec+tp.tv_usec/1e6; }
 	public:
+		static Real getClock(){ timeval tp; gettimeofday(&tp,NULL); return tp.tv_sec+tp.tv_usec/1e6; }
 		Real virtPeriod, virtLast, realPeriod, realLast; long iterPeriod,iterLast,nDo,nDone;
 		PeriodicEngine(): virtPeriod(0),virtLast(0),realPeriod(0),realLast(0),iterPeriod(0),iterLast(0),nDo(-1),nDone(0) { realLast=getClock(); }
 		virtual bool isActivated(){

Modified: trunk/pkg/common/Engine/StandAloneEngine/ResetPositionEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/ResetPositionEngine.cpp	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/pkg/common/Engine/StandAloneEngine/ResetPositionEngine.cpp	2008-12-01 19:47:18 UTC (rev 1587)
@@ -8,7 +8,7 @@
 
 #include"ResetPositionEngine.hpp"
 #include<yade/core/MetaBody.hpp>
-#include<yade/pkg-common/ParticleParameters.hpp>#include<yade/pkg-common/ParticleParameters.hpp>
+#include<yade/pkg-common/ParticleParameters.hpp>
 #include<boost/foreach.hpp>
 
 CREATE_LOGGER(ResetPositionEngine);

Modified: trunk/pkg/common/SConscript
===================================================================
--- trunk/pkg/common/SConscript	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/pkg/common/SConscript	2008-12-01 19:47:18 UTC (rev 1587)
@@ -55,8 +55,8 @@
 		LIBS=env['LIBS']+['Force','ParticleParameters']),
 	env.SharedLibrary('GravityEngines',['Engine/DeusExMachina/GravityEngines.cpp'],
 		LIBS=env['LIBS']+['Force','ParticleParameters']),
-	env.SharedLibrary('JumpChangeSe3',['Engine/DeusExMachina/JumpChangeSe3.cpp'],
-		LIBS=env['LIBS']+['RigidBodyParameters']),
+	env.SharedLibrary('JumpChangeSe3',['Engine/DeusExMachina/JumpChangeSe3.cpp'], LIBS=env['LIBS']+['RigidBodyParameters']),
+	env.SharedLibrary('Se3Interpolator',['Engine/DeusExMachina/Se3Interpolator.cpp'], LIBS=env['LIBS']+['RigidBodyParameters']),
 	env.SharedLibrary('RotationEngine',['Engine/DeusExMachina/RotationEngine.cpp'],
 		LIBS=env['LIBS']+['ParticleParameters','RigidBodyParameters']),
 	env.SharedLibrary('TranslationEngine',['Engine/DeusExMachina/TranslationEngine.cpp'],

Added: trunk/scripts/test/Se3Interpolator.py
===================================================================
--- trunk/scripts/test/Se3Interpolator.py	2008-11-27 08:32:53 UTC (rev 1586)
+++ trunk/scripts/test/Se3Interpolator.py	2008-12-01 19:47:18 UTC (rev 1587)
@@ -0,0 +1,14 @@
+# default parameters or from table
+# encoding: utf-8
+O.bodies.append([
+	utils.sphere([0,0,0],radius=.5),
+])
+
+O.engines=[
+	DeusExMachina('Se3Interpolator',{'subscribedBodies':[0],'goal':[10,10,0,0,0,1,pi],'startIter':10,'goalIter':1010,'goalHook':'print "Finished moving the thing!"; O.pause()'}),
+]
+O.dt=1e-6
+print 'Initial se3:',O.bodies[0].phys['se3']
+#O.saveTmp('init'); O.run(); O.wait();
+#print 'Final   se3:',O.bodies[0].phys['se3']
+#quit()