yade-dev team mailing list archive

[eudoxos at BerliOS <eudoxos@xxxxxxxxxxxxxxx>]

Author: eudoxos
Date: 2009-06-15 17:39:57 +0200 (Mon, 15 Jun 2009)
New Revision: 1798

Modified:
   trunk/gui/py/linterpolation.py
   trunk/gui/qt3/SimulationController.cpp
   trunk/gui/qt3/SimulationController.hpp
   trunk/lib/SConscript
   trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
   trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
   trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
Log:
1. Fix (hopefully) clump damping
2. Remove unused code from InsertionSortCollider
3. Make SimulationController NOT change timestep (rounding) by just opening it (!!)
4. Fix paths in PythonUI and lib/SConscript for $PREFIX/yade$SUFFIX/py/yade instead of $PREFIX/yade$SUFFIX/gui/yade



Modified: trunk/gui/py/linterpolation.py
===================================================================
--- trunk/gui/py/linterpolation.py	2009-06-15 13:43:29 UTC (rev 1797)
+++ trunk/gui/py/linterpolation.py	2009-06-15 15:39:57 UTC (rev 1798)
@@ -73,6 +73,8 @@
 	periodIntegral=integral(x,y)
 	numPeriods=floor(integralValue/periodIntegral)
 	xFrac=xFractionalFromIntegral(integralValue-numPeriods*periodIntegral,x,y)
+	#print '### wanted _%g; period=%g; periodIntegral=_%g (numPeriods=%g); rests _%g (xFrac=%g)'%(integralValue,period,periodIntegral,numPeriods,integralValue-numPeriods*periodIntegral,xFrac)
+	#print '### returning %g*%g+%g=%g'%(period,numPeriods,xFrac,period*numPeriods+xFrac)
 	return period*numPeriods+xFrac
 
 def sanitizeInterpolation(x,y,x0,x1):

Modified: trunk/gui/qt3/SimulationController.cpp
===================================================================
--- trunk/gui/qt3/SimulationController.cpp	2009-06-15 13:43:29 UTC (rev 1797)
+++ trunk/gui/qt3/SimulationController.cpp	2009-06-15 15:39:57 UTC (rev 1798)
@@ -86,6 +86,12 @@
 		loadSimulationFromFileName(Omega::instance().getSimulationFileName());
 	}
 	else{ LOG_DEBUG("Not loading simulation in ctor"); }
+
+	int mantissa, exponent;
+	dtIntegerMantissaExponent(mantissa,exponent);
+	sb10PowerSecond->setValue(exponent);
+	sbSecond->setValue(mantissa);
+
 	// run timer ANY TIME (simulation may be started asynchronously)
 	updateTimerId=startTimer(refreshTime);
 
@@ -352,31 +358,27 @@
 	switch (i)
 	{
 		case 0 : {//Use timeStepper
-			//changeSkipTimeStepper = true;
-			//skipTimeStepper = false;
 			wasUsingTimeStepper=true;
-			//boost::mutex::scoped_lock lock(timeMutex);
 			Omega::instance().skipTimeStepper(false);
 			break;
 			}
 		case 1 : // Try RealTime -- deprecated
-			//changeSkipTimeStepper = true;
-			//skipTimeStepper = true;
-			//wasUsingTimeStepper=false;
 			throw logic_error("RealTime timestep is deprecated and you couldn't click on it!");
 			break;
 		case 2 : {// use fixed time Step
-			//changeSkipTimeStepper = true;
-			//skipTimeStepper = true;
-			//boost::mutex::scoped_lock lock(timeMutex);
 			changeTimeStep = true;
 			wasUsingTimeStepper=false;
 			Omega::instance().skipTimeStepper(true);
 			if(sbSecond->value()==0){ sbSecond->setValue(9); sb10PowerSecond->setValue(sb10PowerSecond->value()-1); }
 			if(sbSecond->value()==10){ sbSecond->setValue(1); sb10PowerSecond->setValue(sb10PowerSecond->value()+1); }
-			Real second = (Real)(sbSecond->value());
-			Real powerSecond = (Real)(sb10PowerSecond->value());
-			Omega::instance().setTimeStep(second*Mathr::Pow(10,powerSecond));
+			int second=(sbSecond->value()), powerSecond = (sb10PowerSecond->value());
+			int exp10,mantissa; dtIntegerMantissaExponent(mantissa,exp10);
+			// only change timestep if the current timestep would have different representation in this integral thing
+			// important so that merely opening the controller doesn't round the existing timestep
+			if((mantissa!=second) || (exp10!=powerSecond)){
+				LOG_DEBUG("Change timestep: current "<<mantissa<<"^"<<exp10<<"; gui "<<second<<"^"<<powerSecond);
+				Omega::instance().setTimeStep((Real)second*Mathr::Pow(10.,(Real)powerSecond));
+			}
 			break;
 		}
 		default: break;
@@ -512,10 +514,13 @@
 	if(rbFixed->isChecked()==usesTimeStepper){ LOG_DEBUG("Checking rbFixed"); rbFixed->setChecked(!usesTimeStepper); }
 	if(rbTimeStepper->isChecked()!=usesTimeStepper){ LOG_DEBUG("Checking rbTimeStepper"); rbTimeStepper->setChecked(usesTimeStepper); }
 
-	Real dt=Omega::instance().getTimeStep();
-	int exp10=floor(log10(dt));
-	sb10PowerSecond->setValue(exp10);
-	sbSecond->setValue((int)(.1+dt/(pow((float)10,exp10)))); // .1: rounding issues
+	int exp10,mantissa; dtIntegerMantissaExponent(mantissa,exp10);
+	sb10PowerSecond->setValue(exp10); sbSecond->setValue(mantissa);
 
 }
 
+void SimulationController::dtIntegerMantissaExponent(int& mantissa, int& exponent){
+	Real dt=Omega::instance().getTimeStep();
+	exponent=floor(log10(dt));
+	mantissa=((int)(.1+dt/(pow((float)10,exponent)))); // .1: rounding issues
+}

Modified: trunk/gui/qt3/SimulationController.hpp
===================================================================
--- trunk/gui/qt3/SimulationController.hpp	2009-06-15 13:43:29 UTC (rev 1797)
+++ trunk/gui/qt3/SimulationController.hpp	2009-06-15 15:39:57 UTC (rev 1798)
@@ -52,6 +52,8 @@
 		QFrame * scrollViewFrame;
 		QVBoxLayout* scrollViewLayout;
 		void addNewView();
+		void dtIntegerMantissaExponent(int& mantissa, int& exponent);
+
 	
 	public : 
 		void loadSimulationFromFileName(const std::string& fileName,bool center=true);

Modified: trunk/lib/SConscript
===================================================================
--- trunk/lib/SConscript	2009-06-15 13:43:29 UTC (rev 1797)
+++ trunk/lib/SConscript	2009-06-15 15:39:57 UTC (rev 1798)
@@ -2,7 +2,7 @@
 Import('*')
 
 if 'EMBED_PYTHON' in env['CPPDEFINES']:
-	env.Install('$PREFIX/lib/yade$SUFFIX/gui/yade',[
+	env.Install('$PREFIX/lib/yade$SUFFIX/py/yade',[
 		env.SharedLibrary('WeightedAverage2d',['smoothing/WeightedAverage2d.cpp'],SHLIBPREFIX='',
 		LIBS=env['LIBS']+['yade-base'])
 	])

Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp	2009-06-15 13:43:29 UTC (rev 1797)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp	2009-06-15 15:39:57 UTC (rev 1798)
@@ -111,20 +111,6 @@
 
 	// process interactions that the constitutive law asked to be erased
 	interactions->erasePending(*this);
-	#if 0
-	FOREACH(const InteractionContainer::bodyIdPair& p, interactions->pendingErase){
-		// remove those that do not overlap spatially anymore
-		if(!spatialOverlap(p[0],p[1])){ interactions->erase(p[0],p[1]); LOG_TRACE("Deleted interaction #"<<p[0]<<"+#"<<p[1]); }
-		else
-		{
-			const shared_ptr<Interaction>& I=interactions->find(p[0],p[1]);
-			if(!I){ LOG_FATAL("Requested deletion of a non-existent interaction #"<<p[0]<<"+#"<<p[1]<<"?!"); throw; }
-			I->reset();
-		}
-	}
-	interactions->pendingErase.clear();
-	#endif
-	
 
 	// sort
 		if(!doInitSort && !sortThenCollide){
@@ -134,6 +120,7 @@
 		else {
 			if(doInitSort){
 				// the initial sort is in independent in 3 dimensions, may be run in parallel
+				// it seems that there is no time gain running this in parallel, though
 				#pragma omp parallel sections
 				{
 					#pragma omp section
@@ -159,7 +146,7 @@
 				// go up until we meet the upper bound
 				for(size_t j=i+1; V[j].id!=iid; j++){
 					const body_id_t& jid=V[j].id;
-					/// FIXME: not sure why this doesn't work. If this condition is commented out, we have exact same interactions as from SpatialQuickSort. Otherwise some interactions are missing!
+					/// Not sure why this doesn't work. If this condition is commented out, we have exact same interactions as from SpatialQuickSort. Otherwise some interactions are missing!
 					// skip bodies with smaller (arbitrary, could be greater as well) id, since they will detect us when their turn comes
 					// if(jid<iid) { /* LOG_TRACE("Skip #"<<V[j].id<<(V[j].flags.isMin?"(min)":"(max)")<<" with "<<iid<<" (smaller id)"); */ continue; }
 					/* abuse the same function here; since it does spatial overlap check first, it is OK to use it */

Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp	2009-06-15 13:43:29 UTC (rev 1797)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp	2009-06-15 15:39:57 UTC (rev 1798)
@@ -26,20 +26,27 @@
 	damping = 0.2;
 }
 
+void NewtonsDampedLaw::cundallDamp(const Real& dt, const Vector3r& f, const Vector3r& velocity, Vector3r& acceleration, const Vector3r& m, const Vector3r& angularVelocity, Vector3r& angularAcceleration){
+	for(int i=0; i<3; i++){
+		angularAcceleration[i]*= 1 - damping*Mathr::Sign ( m[i]*(angularVelocity[i] + (Real) 0.5 *dt*angularAcceleration[i]) );
+		acceleration       [i]*= 1 - damping*Mathr::Sign ( f[i]*(velocity       [i] + (Real) 0.5 *dt*acceleration       [i]) );
+	}
+}
+
 void NewtonsDampedLaw::applyCondition ( MetaBody * ncb )
 {
 	ncb->bex.sync();
+	Real dt=Omega::instance().getTimeStep();
+
 	FOREACH(const shared_ptr<Body>& b, *ncb->bodies){
 		// clump members are non-dynamic; they skip the rest of loop once their forces are properly taken into account, however
 		if (!b->isDynamic && !b->isClumpMember()) continue;
 		
 		RigidBodyParameters* rb = YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
-		body_id_t id = b->getId();
+		const body_id_t& id=b->getId();
 		const Vector3r& m=ncb->bex.getTorque(id);
 		const Vector3r& f=ncb->bex.getForce(id);
 
-		Real dt = Omega::instance().getTimeStep();
-
 		//Newtons mometum law :
 		if (b->isStandalone()){
 			rb->angularAcceleration=diagDiv(m,rb->inertia);
@@ -53,27 +60,23 @@
 		else {
 			assert(b->isClumpMember());
 			assert(b->clumpId>b->id);
-			const shared_ptr<Body>& clump ( Body::byId ( b->clumpId ) );
+			const shared_ptr<Body>& clump=Body::byId(b->clumpId,ncb);
 			RigidBodyParameters* clumpRBP=YADE_CAST<RigidBodyParameters*> ( clump->physicalParameters.get() );
-			/* angularAcceleration is reset by ClumpMemberMover engine */
-			clumpRBP->angularAcceleration+=diagDiv ( m,clumpRBP->inertia );
-			clumpRBP->acceleration+=f/clumpRBP->mass;
-			clumpRBP->angularAcceleration+=diagDiv ( ( rb->se3.position-clumpRBP->se3.position ).Cross ( f ),clumpRBP->inertia ); //acceleration from torque generated by the force WRT particle centroid on the clump centroid
+			Vector3r diffClumpAccel=f/clumpRBP->mass;
+			// angular acceleration from: normal torque + torque generated by the force WRT particle centroid on the clump centroid
+			Vector3r diffClumpAngularAccel=diagDiv(m,clumpRBP->inertia)+diagDiv((rb->se3.position-clumpRBP->se3.position).Cross(f),clumpRBP->inertia); 
+			// damp increment of accels on the clump, using velocities of the clump MEMBER
+			cundallDamp(dt,f,rb->velocity,diffClumpAccel,m,rb->angularVelocity,diffClumpAngularAccel);
+			// clumpRBP->{acceleration,angularAcceleration} are reset byt Clump::moveMembers, it is ok to just increment here
+			clumpRBP->acceleration+=diffClumpAccel;
+			clumpRBP->angularAcceleration+=diffClumpAngularAccel;
 			continue;
 		}
 
+		assert(!b->isClumpMember());
+		// damping: applied to non-clumps only, as clumps members were already damped above
+		if(!b->isClump()) cundallDamp(dt,f,rb->velocity,rb->acceleration,m,rb->angularVelocity,rb->angularAcceleration);
 
-		// damping: applied to both clumps and non-clumps
-		// FIXME: clumps with just 1 body are not dumped? See sphere jumping forever in scripts/test/clump.py
-		//   other clumps seem to be damped in some strange way (again, see the same script).
-		//   The damping also looks different from ClumpTestGen, which uses separate engines instead of NewtonsDampedLaw
-		for ( int i=0; i<3; ++i ) {
-			//Damping moments
-			rb->angularAcceleration[i] *= 1 - damping*Mathr::Sign ( m[i]*(rb->angularVelocity[i] + (Real) 0.5 *dt*rb->angularAcceleration[i]) );
-			//Damping force :
-			rb->acceleration[i] *= 1 - damping*Mathr::Sign ( f[i]*(rb->velocity[i] + (Real) 0.5 *dt*rb->acceleration[i]) );
-		}
-
 		// blocking DOFs
 		if(rb->blockedDOFs==0){ /* same as: rb->blockedDOFs==PhysicalParameters::DOF_NONE */
 			rb->angularVelocity=rb->angularVelocity+dt*rb->angularAcceleration;

Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp	2009-06-15 13:43:29 UTC (rev 1797)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp	2009-06-15 15:39:57 UTC (rev 1798)
@@ -35,6 +35,7 @@
 
 class NewtonsDampedLaw : public DeusExMachina
 {
+	inline void cundallDamp(const Real& dt, const Vector3r& f, const Vector3r& velocity, Vector3r& acceleration, const Vector3r& m, const Vector3r& angularVelocity, Vector3r& angularAcceleration);
 	public :
 		///damping coefficient for Cundall's non viscous damping
 		Real damping;