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;
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
