Author: eudoxos
Date: 2009-07-13 21:31:44 +0200 (Mon, 13 Jul 2009)
New Revision: 1862
Added:
trunk/pkg/common/DataClass/VelocityBins.cpp
trunk/pkg/common/DataClass/VelocityBins.hpp
Removed:
trunk/pkg/common/DataClass/State/
Modified:
trunk/gui/qt3/QtGUI.cpp
trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.cpp
trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.hpp
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
trunk/pkg/common/SConscript
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
trunk/pkg/dem/SConscript
trunk/scripts/test/collider-stride-triax.py
Log:
1. Add velocity binning for collision optimization: requires cooperation of
NewtonsDampedLaw, BoundingVolumeMetaEngine and InsertionSortCollider.
2. Remove time-based sweeping from BoundingVOlumeMetaEngine.
3. QtGUI doesn't show annoying warning when run from python session anymore.
Modified: trunk/gui/qt3/QtGUI.cpp
===================================================================
--- trunk/gui/qt3/QtGUI.cpp 2009-07-13 15:30:36 UTC (rev 1861)
+++ trunk/gui/qt3/QtGUI.cpp 2009-07-13 19:31:44 UTC (rev 1862)
@@ -84,9 +84,14 @@
void QtGUI::runNaked(){
if(!app){ // no app existing yet
+ if(getenv("DISPLAY")==0){
+ LOG_ERROR("$DISPLAY environment var not set, not
starting qt3 gui.");
+ return;
+ };
LOG_INFO("Creating QApplication");
XInitThreads();
- app=new QApplication(0,NULL);
+ int argc=0;
+ app=new QApplication(argc,(char**)NULL);
if(!YadeQtMainWindow::self){
YadeQtMainWindow::guiMayDisappear=true;
mainWindow=new YadeQtMainWindow();
Added: trunk/pkg/common/DataClass/VelocityBins.cpp
===================================================================
--- trunk/pkg/common/DataClass/VelocityBins.cpp 2009-07-13 15:30:36 UTC (rev
1861)
+++ trunk/pkg/common/DataClass/VelocityBins.cpp 2009-07-13 19:31:44 UTC (rev
1862)
@@ -0,0 +1,97 @@
+// 2009 © Václav Šmilauer <[email protected]>
+
+
+#include"VelocityBins.hpp"
+#include<yade/core/MetaBody.hpp>
+#include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<boost/foreach.hpp>
+#ifndef FOREACH
+# define FOREACH BOOST_FOREACH
+#endif
+CREATE_LOGGER(VelocityBins);
+
+bool VelocityBins::incrementDists_shouldCollide(Real dt){
+ int i=0;
+ FOREACH(Bin& bin, bins){
+ bin.currDistSq+=dt*dt*bin.currMaxVelSq; i++;
+ if(bin.currDistSq>pow(bin.maxDist,2)){
+ LOG_TRACE("Collide: bin"<<i<<": max dist
"<<bin.maxDist<<", current "<<sqrt(bin.currDistSq));
+ return true;
+ }
+ }
+ return false;
+}
+
+void VelocityBins::setBins(MetaBody* rootBody, Real currMaxVelSq, Real
refSweepLength){
+ // initialization
+ // sanity checks
+ if(nBins<1 || nBins>100){LOG_FATAL("Number of bins must be >=1
and <=100"); abort(); }
+ if(binOverlap>=1 || binOverlap<=0){ LOG_ERROR("binOverlap set
to 0.8 (was "<<binOverlap<<", not in range (0…1) )"); binOverlap=0.8;}
+ // number of bins changed
+ if(nBins!=bins.size()){
+ LOG_INFO("New number of bins: "<<nBins);
+ bins.resize(nBins);
+ }
+ // number of bodies changed
+ if(bodyBins.size()!=rootBody->bodies->size())
bodyBins.resize(rootBody->bodies->size(),-1);
+ // set the new overall refMaxVelSq
+ if(refMaxVelSq<0){ refMaxVelSq=currMaxVelSq; /* first time */}
+ else {
+ // there should be some maximum speed change parameter,
so that bins do not change their limits (and therefore bodies, also!) too
often, depending on 1 particle going crazy
+
refMaxVelSq=min(max(refMaxVelSq*pow(1-maxRefRelStep,2),currMaxVelSq),refMaxVelSq*pow(1+maxRefRelStep,2));
+ if(refMaxVelSq==0) refMaxVelSq=currMaxVelSq;
+ }
+ LOG_TRACE("new refMaxVel: "<<sqrt(refMaxVelSq));
+ // compute new minima/maxima for all bins, reset distance
counters and real max velocity
+ for(size_t i=0; i<nBins; i++){
+ Bin& bin=bins[i];
+ // 0th bin (fastest) has maximum the current maximum;
slowest bin has minimum 0.
+ bin.binMaxVelSq=(i==0 ? currMaxVelSq :
refMaxVelSq/pow(binCoeff*binCoeff,i));
+ bin.binMinVelSq=(i==nBins-1 ? 0. :
refMaxVelSq/pow(binCoeff*binCoeff,i+1));
+ bin.maxDist=(i==0 ?
+ (refMaxVelSq==0 ? 0:
sqrt(currMaxVelSq/refMaxVelSq)*abs(refSweepLength)) :
+ abs(refSweepLength)/pow(binCoeff,i)
+ );
+ bin.currDistSq=0; bin.currMaxVelSq=0; bin.nBodies=0;
+ }
+ long moveFaster=0, moveSlower=0;
+ FOREACH(const shared_ptr<Body>& b, *rootBody->bodies){
+ RigidBodyParameters*
rbp=YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
+ Real velSq=VelocityBins::getBodyVelSq(rbp);
+ binNo_t newBin=-1, oldBin=bodyBins[b->getId()];
+ // we could compute logarithm and round it here, but perhaps
this is faster
+ for(size_t i=0; i<nBins; i++){
+ // 1. do not allow to move up (slow down) more than one
bin at time (reasonable?)
+ // 2. for the current bin, put the lower margin lower
by binOverlap^2, to avoid too many oscillations
+ if((oldBin>=0 && oldBin==(binNo_t)(i-1)) ||
velSq>=bins[i].binMinVelSq*(oldBin==(binNo_t)i ? pow(binOverlap,2) : 1.) ){
newBin=(binNo_t)i; break;}
+ }
+ if(newBin<0){ LOG_FATAL("Body #"<<b->getId()<<",
velSq="<<velSq<<" not put in any bin!?"); abort(); }
+ if(oldBin>=0) { if(newBin>oldBin) moveSlower++; else
if(oldBin>newBin) moveFaster++; }
+ bodyBins[b->getId()]=newBin;
+ // LOG_TRACE("#"<<b->getId()<<": vel="<<sqrt(velSq)<<", bin
"<<(int)newBin);
+ bins[newBin].nBodies+=1;
+ }
+ #ifdef LOG4CXX
+ // if debugging output
+ if(logger->isDebugEnabled() &&
(rootBody->currentIteration-histLast>=histInterval || histLast<0)){
+ histLast=rootBody->currentIteration;
+ LOG_INFO(bodyBins.size()<<" bodies (moves:
"<<moveFaster<<" faster, "<<moveSlower<<" slower),
refMaxVel="<<sqrt(refMaxVelSq));
+ for(size_t i=0; i<nBins; i++){
+ int
nChars=int(80*(bins[i].nBodies/Real(bodyBins.size())));
+ string l; for(int j=0; j<nChars; j++) l+="#";
for(int j=nChars; j<80; j++) l+=" ";
+ cerr<<"\t"<<i<<": |"<<l<<"|
("<<bins[i].nBodies<<")
"<<sqrt(bins[i].binMinVelSq)<<"…"<<sqrt(bins[i].binMaxVelSq)<<endl;
+ }
+ }
+ #endif
+}
+
+/* non-parallel implementations */
+void VelocityBins::binVelSqInitialize(){ FOREACH(Bin& bin, bins)
bin.currMaxVelSq=0; }
+void VelocityBins::binVelSqUse(body_id_t id, Real velSq){
+ #ifdef YADE_OPENMP
+ if(omp_get_thread_num()>0) {LOG_FATAL("VelocityBins do not
support multiple openMP threads yet!"); abort(); }
+ #endif
+ Real& maxVelSq(bins[bodyBins[id]].currMaxVelSq);
+ maxVelSq=max(maxVelSq,velSq);
+}
+void VelocityBins::binVelSqFinalize(){}
Added: trunk/pkg/common/DataClass/VelocityBins.hpp
===================================================================
--- trunk/pkg/common/DataClass/VelocityBins.hpp 2009-07-13 15:30:36 UTC (rev
1861)
+++ trunk/pkg/common/DataClass/VelocityBins.hpp 2009-07-13 19:31:44 UTC (rev
1862)
@@ -0,0 +1,76 @@
+// 2009 © Václav Šmilauer <[email protected]>
+#pragma once
+
+#include<yade/core/Interaction.hpp> // for body_id_t
+#include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<vector>
+
+class MetaBody;
+
+/* Class for putting bodies to velcoity bins, for optimization
+of collision detection.
+Each bin is characterized by its minimum/maximum velocity.
+See
http://yade.wikia.com/wiki/Insertion_Sort_Collider_Stride#Enhancement_ideas:_velocity_bins
for brief design overview.
+*/
+class VelocityBins{
+ public:
+ VelocityBins(int _nBins, Real _refMaxVelSq, Real _binCoeff=10, Real
_binOverlap=0.8): refMaxVelSq(_refMaxVelSq), binCoeff(_binCoeff),
binOverlap(_binOverlap), maxRefRelStep(0.05), nBins(_nBins), histInterval(200),
histLast(-1){};
+ typedef signed char binNo_t;
+ struct Bin{
+ Bin(): binMinVelSq(-1), binMaxVelSq(-1), maxDist(0),
currDistSq(0), currMaxVelSq(0), nBodies(0){}
+ // limits for bin memebrship
+ Real binMinVelSq, binMaxVelSq;
+ // maximum distance that body in this bin can travel before it
goes out of its swept bbox
+ Real maxDist;
+ // distance so far traveled by the fastest body in this bin
(since last setBins)
+ Real currDistSq;
+ // maximum velSq over all bodies in this bin
+ Real currMaxVelSq;
+ // number of bodies in this bin (for informational purposes
only)
+ long nBodies;
+ };
+ // bins themselves (their number is nBins)
+ vector<Bin> bins;
+ // each body has its bin number stored here
+ vector<binNo_t> bodyBins;
+ // reference overall maximum velocity (may be different from
bins[0].maxVelSq)
+ Real refMaxVelSq;
+ // factor by which maximum velocity in each bin is smaller than in the
higer one
+ Real binCoeff;
+ // relative overlap beween bins; body will not be moved from faster bin
until its velocity is min*binOverlap; must be <=1
+ Real binOverlap;
+ // maximum relative change of reference max velocity per invocation
+ Real maxRefRelStep;
+ // number of bins; must be >=1 and <=100 (artificial upper limit)
+ size_t nBins;
+ // how often to show histogram, if LOG_DEBUG is enabled;
+ long histInterval, histLast;
+ // Assign bins to all bodies
+ void setBins(MetaBody*, Real currMaxVelSq, Real refSweepLength);
+
+ // Increment maximum per-bin distances and tell whether some bodies may
be already getting out of the swept bbox (in that case, we need to
recompute bounding volumes and run the collider)
+ bool incrementDists_shouldCollide(Real dt);
+
+ /* NOTE: following 3 functions are separated because of multi-threaded
operation of NewtonsDampedLaw
+ in that case, every thread must have its own per-bin maximum and
binVelSqFinalize will assign the
+ overall maxima to the bins.
+
+ The non-parallel implementation simply resets all Bin::realMaxVelSq,
maximize if for each
+ body in that bin and binVelSqFinalize() will not do nothing.
+ */
+ // reset per-bin max velocities
+ void binVelSqInitialize();
+ // use body max velocity -- called for every body at every step (from
NewtonsDampedLaw, normally)
+ void binVelSqUse(body_id_t id, Real velSq);
+ // actually assign max velocities to their respective bins
+ void binVelSqFinalize();
+
+ // get velSq for given body; this should be called from
NewtonsDampedLaw as well,
+ // to ensure that the same formulas are used (once we have
angularVelocity + AABB span,
+ // for instance
+ static Real getBodyVelSq(RigidBodyParameters* rbp){
+ return rbp->velocity.SquaredLength();
+ }
+
+ DECLARE_LOGGER;
+};
Modified: trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.cpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.cpp
2009-07-13 15:30:36 UTC (rev 1861)
+++ trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.cpp
2009-07-13 19:31:44 UTC (rev 1862)
@@ -14,6 +14,7 @@
#include<yade/core/MetaBody.hpp>
#include<yade/pkg-common/AABB.hpp>
#include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<yade/pkg-common/VelocityBins.hpp>
CREATE_LOGGER(BoundingVolumeMetaEngine);
@@ -21,6 +22,8 @@
{
shared_ptr<BodyContainer>& bodies = ncb->bodies;
const long numBodies=(long)bodies->size();
+ bool haveBins=(bool)velocityBins;
+ if(sweepDist!=0 && haveBins){ LOG_FATAL("Only one of sweepDist or
velocityBins can used!"); abort(); }
//#pragma omp parallel for
for(int id=0; id<numBodies; id++){
if(!bodies->exists(id)) continue; // don't delete this check -
Janek
@@ -36,28 +39,15 @@
#else
operator()(ig,b->boundingVolume,b->physicalParameters->se3,b.get());
#endif
- if(sweepTime>0){
- /* compute se3 the body would have after given time
(linear interpolation), and use the maximum bounding volume then.
- It is a pessimisation that allows us to not update the
aabb during some period of time. */
- Vector3r min0=b->boundingVolume->min,
max0=b->boundingVolume->max;
- const Se3r& se3=b->physicalParameters->se3;
- const RigidBodyParameters*
rbp=YADE_CAST<RigidBodyParameters*>(b->physicalParameters.get());
- /* For position change, Δr=vΔt+½aΔt²; similar for
orientation. */
- Se3r sweptSe3;
sweptSe3.position=se3.position+(rbp->velocity+.5*rbp->acceleration*sweepTime)*sweepTime;
- Vector3r
axis=rbp->angularVelocity+.5*rbp->angularAcceleration*sweepTime; Real
angle=axis.Normalize(); Quaternionr q(axis,angle*sweepTime);
sweptSe3.orientation=q*se3.orientation;
- #ifdef BV_FUNCTOR_CACHE
-
ig->boundFunctor->go(ig,b->boundingVolume,sweptSe3,b.get());
- #else
-
operator()(ig,b->boundingVolume,sweptSe3,b.get());
- #endif
- AABB* aabb=YADE_CAST<AABB*>(b->boundingVolume.get());
// to update halfSize and center as well
- aabb->min=componentMinVector(min0,aabb->min);
aabb->max=componentMaxVector(max0,aabb->max);
- aabb->center=.5*(aabb->min+aabb->max);
aabb->halfSize=.5*(aabb->max-aabb->min);
+ if(sweepDist>0){
+ AABB* aabb=YADE_CAST<AABB*>(b->boundingVolume.get());
+ aabb->halfSize+=Vector3r(sweepDist,sweepDist,sweepDist);
+ aabb->min=aabb->center-aabb->halfSize;
aabb->max=aabb->center+aabb->halfSize;
}
- if(sweepDist!=0){
+ if(haveBins){
AABB* aabb=YADE_CAST<AABB*>(b->boundingVolume.get());
- if(sweepDist<0)
aabb->halfSize=-sweepDist*aabb->halfSize; // relative scaling by abs(sweepDist)
- else
aabb->halfSize+=Vector3r(sweepDist,sweepDist,sweepDist); // absolute distance
increment
+ Real
sweep=velocityBins->bins[velocityBins->bodyBins[b->getId()]].maxDist;
+ aabb->halfSize+=Vector3r(sweep,sweep,sweep);
aabb->min=aabb->center-aabb->halfSize;
aabb->max=aabb->center+aabb->halfSize;
}
}
Modified: trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.hpp
===================================================================
--- trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.hpp
2009-07-13 15:30:36 UTC (rev 1861)
+++ trunk/pkg/common/Engine/MetaEngine/BoundingVolumeMetaEngine.hpp
2009-07-13 19:31:44 UTC (rev 1862)
@@ -16,6 +16,8 @@
#include<yade/core/BoundingVolume.hpp>
#include<yade/core/Body.hpp>
+class VelocityBins;
+
class BoundingVolumeMetaEngine : public MetaEngine2D
<
InteractingGeometry,
// base classe for dispatch
@@ -34,19 +36,16 @@
public :
// selectively turn off this engine (from within the collider,
for instance)
bool activated;
- //! if < 0, will set bounding volumes that the bodies (most
likely) will not go over within given time interval; otherwise, do not sweep
- Real sweepTime;
- //! bounding box will be enlarged by this amount in all 3
dimensions;
- /// if negative, it is relative size (e.g. -2 will make the
bounding box 2 × bigger)
- /// if positive, it is absolute distance; its will be added in
all 6 directions (i.e. dimensions will be enlarged by 2 × this value)
+ //! bounding box will be enlarged by this amount in all 3
dimensions; must be non-negative;
Real sweepDist;
- BoundingVolumeMetaEngine(): activated(true), sweepTime(-1),
sweepDist(0) {}
+ BoundingVolumeMetaEngine(): activated(true), sweepDist(0) {}
virtual void action(MetaBody*);
virtual bool isActivated(MetaBody*){ return activated; }
+ shared_ptr<VelocityBins> velocityBins;
DECLARE_LOGGER;
REGISTER_CLASS_NAME(BoundingVolumeMetaEngine);
REGISTER_BASE_CLASS_NAME(MetaEngine2D);
- REGISTER_ATTRIBUTES(MetaEngine,(activated)(sweepTime)(sweepDist));
+ REGISTER_ATTRIBUTES(MetaEngine,(activated)(sweepDist));
};
REGISTER_SERIALIZABLE(BoundingVolumeMetaEngine);
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
2009-07-13 15:30:36 UTC (rev 1861)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.cpp
2009-07-13 19:31:44 UTC (rev 1862)
@@ -5,6 +5,7 @@
#include<yade/core/Interaction.hpp>
#include<yade/core/InteractionContainer.hpp>
#include<yade/pkg-common/BoundingVolumeMetaEngine.hpp>
+#include<yade/pkg-common/VelocityBins.hpp>
#include<yade/pkg-dem/NewtonsDampedLaw.hpp>
#include<algorithm>
@@ -65,12 +66,15 @@
bool InsertionSortCollider::isActivated(MetaBody* rb){
// activated if number of bodies changes (hence need to refresh
collision information)
// or the time of scheduled run already came, or we were never
scheduled yet
- if(sweepLength<=0) return true;
- if(!newton || newton->maxVelocitySq<0 || sweepVelocity<0)
return true; // we wouldn't be able to find the max velocity, or it has not
been computed yet, or there were no data available last time
+ if(!strideActive) return true;
+ if(!newton || (nBins>=1 && !newton->velocityBins)) return true;
+ if(nBins>=1 &&
newton->velocityBins->incrementDists_shouldCollide(rb->dt)) return true;
+ if(nBins<=0){
+ if(fastestBodyMaxDist<0){fastestBodyMaxDist=0; return
true;}
+ fastestBodyMaxDist+=sqrt(newton->maxVelocitySq)*rb->dt;
+ if(fastestBodyMaxDist>=sweepLength) return true;
+ }
if(XX.size()!=2*rb->bodies->size()) return true;
- if(fastestBodyMaxDist<0){fastestBodyMaxDist=0; return true;}
- fastestBodyMaxDist+=sqrt(newton->maxVelocitySq)*rb->dt;
- if(fastestBodyMaxDist>=sweepLength) return true;
// we wouldn't run in this step; in that case, just delete
pending interactions
// this is done in ::action normally, but it would make the
call counters not reflect the stride
rb->interactions->erasePending(*this);
@@ -128,21 +132,37 @@
ISC_CHECKPOINT("init");
#ifdef COLLIDE_STRIDED
- if(sweepLength>0){
- if(newton->maxVelocitySq>=0){ // non-negative,
i.e. a really computed value
- // compute new stride value
- assert(sweepFactor>1.);
-
sweepVelocity=sqrt(newton->maxVelocitySq)*sweepFactor;
+ // get us ready for strides
+ if(!strideActive && sweepLength>0){
+ if(newton->maxVelocitySq>=0){ // maxVelocitySq
is a really computed value
+ strideActive=true;
+ if(nBins>=1){
+ if(!newton->velocityBins){
newton->velocityBins=shared_ptr<VelocityBins>(new
VelocityBins(nBins,newton->maxVelocitySq,binCoeff,binOverlap)); }
+
if(!boundDispatcher->velocityBins)
boundDispatcher->velocityBins=newton->velocityBins;
+ }
+ }
+ }
+ if(strideActive){
+ assert(sweepLength>0);
+ if(nBins<=0){
+ // reset bins, in case they were active
but are not anymore
+ if(newton->velocityBins)
newton->velocityBins=shared_ptr<VelocityBins>();
if(boundDispatcher->velocityBins)
boundDispatcher->velocityBins=shared_ptr<VelocityBins>();
+ assert(strideActive);
assert(newton->maxVelocitySq>=0); assert(sweepFactor>1.);
+ Real
sweepVelocity=sqrt(newton->maxVelocitySq)*sweepFactor; int stride=-1;
if(sweepVelocity>0) {
stride=max(1,int((sweepLength/sweepVelocity)/rb->dt));
boundDispatcher->sweepDist=rb->dt*(stride-1)*sweepVelocity;
} else { // no motion
boundDispatcher->sweepDist=0;
// nothing moves, no need to make bboxes larger
}
- } else { /* no valid data yet, run next time
again */ boundDispatcher->sweepDist=0; sweepVelocity=-1; stride=1; }
- LOG_DEBUG(rb->simulationTime<<"s: stride
adapted to "<<stride<<"; sweepVelocity="<<sweepVelocity<<",
maxVelocity="<<sqrt(newton->maxVelocitySq)<<",
sweepDist="<<boundDispatcher->sweepDist);
- fastestBodyMaxDist=0; // reset
- } else { boundDispatcher->sweepTime=-1;
boundDispatcher->sweepDist=0; }
+ LOG_DEBUG(rb->simulationTime<<"s:
stride ≈"<<stride<<"; maxVelocity="<<sqrt(newton->maxVelocitySq)<<",
sweepDist="<<boundDispatcher->sweepDist);
+ fastestBodyMaxDist=0; // reset
+ } else { // nBins>=1
+ assert(newton->velocityBins);
assert(boundDispatcher->velocityBins);
+ // re-bin bodies
+
newton->velocityBins->setBins(rb,newton->maxVelocitySq,sweepLength);
+ }
+ }
boundDispatcher->action(rb);
#endif
Modified: trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
===================================================================
--- trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
2009-07-13 15:30:36 UTC (rev 1861)
+++ trunk/pkg/common/Engine/StandAloneEngine/InsertionSortCollider.hpp
2009-07-13 19:31:44 UTC (rev 1862)
@@ -46,21 +46,22 @@
shared_ptr<BoundingVolumeMetaEngine> boundDispatcher;
// we need this to find out about current maxVelocitySq
shared_ptr<NewtonsDampedLaw> newton;
- // interval at which we will run (informative only, is updated
automatically)
- int stride;
+ // if False, no type of striding is used
+ // if True, then either sweepLength XOR nBins is set
+ bool strideActive;
/// Absolute length that will be added to bounding boxes at
each side; it should be something like 1/5 of typical grain radius
/// this value is used to adapt stride; if too large, stride
will be big, but the ratio of potential-only interactions will be very big,
/// thus slowing down collider & interaction loops
significantly (remember: O(addLength^3))
/// If non-positive, collider runs always, without stride
adaptivity
Real sweepLength;
- //! Current sweeping velocity; computed from
maxVelocity*sweepVelocityFactor
- Real sweepVelocity;
//! Overestimation factor for the sweep velocity; must be >=1.0.
/// Has no influence on sweepLength, only on the computed
stride.
/// Default 1.05
Real sweepFactor;
//! maximum distance that the fastest body could have travelled
since the last run; if >= sweepLength, we could get out of bboxes and will
trigger full run
Real fastestBodyMaxDist;
+ // parameters to be passed to VelocityBins, if nBins>0
+ int nBins; Real binCoeff, binOverlap;
#endif
//! storage for bounds
std::vector<Bound> XX,YY,ZZ;
@@ -90,7 +91,7 @@
InsertionSortCollider():
#ifdef COLLIDE_STRIDED
- stride(0), sweepLength(-1), sweepVelocity(-1),
sweepFactor(1.05), fastestBodyMaxDist(-1),
+ strideActive(false), sweepLength(-1), sweepFactor(1.05),
fastestBodyMaxDist(-1), nBins(0), binCoeff(5), binOverlap(0.8),
#endif
sortAxis(0), sortThenCollide(false){
#ifdef ISC_TIMING
@@ -101,7 +102,7 @@
REGISTER_CLASS_AND_BASE(InsertionSortCollider,Collider);
REGISTER_ATTRIBUTES(Collider,(sortAxis)(sortThenCollide)
#ifdef COLLIDE_STRIDED
-
(stride)(sweepLength)(sweepFactor)(sweepVelocity)(fastestBodyMaxDist)
+
(strideActive)(sweepLength)(sweepFactor)(fastestBodyMaxDist)(nBins)(binCoeff)(binOverlap)
#endif
);
DECLARE_LOGGER;
Modified: trunk/pkg/common/SConscript
===================================================================
--- trunk/pkg/common/SConscript 2009-07-13 15:30:36 UTC (rev 1861)
+++ trunk/pkg/common/SConscript 2009-07-13 19:31:44 UTC (rev 1862)
@@ -10,8 +10,8 @@
env.Install('$PREFIX/lib/yade$SUFFIX/pkg-common',[
- env.SharedLibrary('ParallelEngine',
- ['Engine/ParallelEngine.cpp']),
+ env.SharedLibrary('VelocityBins',['DataClass/VelocityBins.cpp']), # not
a plugin
+ env.SharedLibrary('ParallelEngine',['Engine/ParallelEngine.cpp']),
env.SharedLibrary('EngineUnits',['Engine/MetaEngine/EngineUnits.cpp']),
env.SharedLibrary('AABB',['DataClass/BoundingVolume/AABB.cpp']),
env.SharedLibrary('BoundingSphere',['DataClass/BoundingVolume/BoundingSphere.cpp']),
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-13
15:30:36 UTC (rev 1861)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.cpp 2009-07-13
19:31:44 UTC (rev 1862)
@@ -9,6 +9,7 @@
#include"NewtonsDampedLaw.hpp"
#include<yade/core/MetaBody.hpp>
#include<yade/pkg-common/RigidBodyParameters.hpp>
+#include<yade/pkg-common/VelocityBins.hpp>
#include<yade/lib-base/yadeWm3Extra.hpp>
#include<yade/pkg-dem/Clump.hpp>
@@ -34,6 +35,7 @@
// clumpRBP->{acceleration,angularAcceleration} are reset byt
Clump::moveMembers, it is ok to just increment here
clumpRBP->acceleration+=diffClumpAccel;
clumpRBP->angularAcceleration+=diffClumpAngularAccel;
+ if(haveBins)
velocityBins->binVelSqUse(memberId,VelocityBins::getBodyVelSq(rb));
maxVelocitySq=max(maxVelocitySq,rb->velocity.SquaredLength());
}
@@ -42,6 +44,8 @@
ncb->bex.sync();
Real dt=Omega::instance().getTimeStep();
maxVelocitySq=-1;
+ haveBins=(bool)velocityBins;
+ if(haveBins) velocityBins->binVelSqInitialize();
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
@@ -67,6 +71,7 @@
rb->angularAcceleration+=dAngAccel;
}
+ if(haveBins)
{velocityBins->binVelSqUse(id,VelocityBins::getBodyVelSq(rb));}
maxVelocitySq=max(maxVelocitySq,rb->velocity.SquaredLength());
// blocking DOFs
@@ -94,5 +99,7 @@
if(b->isClump()) static_cast<Clump*>(b.get())->moveMembers();
}
+
+ if(haveBins) velocityBins->binVelSqFinalize();
}
Modified: trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp
===================================================================
--- trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-13
15:30:36 UTC (rev 1861)
+++ trunk/pkg/dem/Engine/DeusExMachina/NewtonsDampedLaw.hpp 2009-07-13
19:31:44 UTC (rev 1862)
@@ -32,18 +32,22 @@
*/
class RigidBodyParameters;
+class VelocityBins;
class NewtonsDampedLaw : public StandAloneEngine{
inline void cundallDamp(const Real& dt, const Vector3r& f, const
Vector3r& velocity, Vector3r& acceleration, const Vector3r& m, const Vector3r&
angularVelocity, Vector3r& angularAcceleration);
void handleClumpMember(MetaBody* ncb, const body_id_t memberId,
RigidBodyParameters* clumpRBP);
-
+ bool haveBins;
public:
///damping coefficient for Cundall's non viscous damping
Real damping;
/// store square of max. velocity, for informative purposes;
computed again at every step
Real maxVelocitySq;
+ /// velocity bins (not used if not created)
+ shared_ptr<VelocityBins> velocityBins;
virtual void action(MetaBody *);
NewtonsDampedLaw(): damping(0.2), maxVelocitySq(-1){}
+
REGISTER_ATTRIBUTES(StandAloneEngine,(damping)(maxVelocitySq));
REGISTER_CLASS_AND_BASE(NewtonsDampedLaw,StandAloneEngine);
};
Modified: trunk/pkg/dem/SConscript
===================================================================
--- trunk/pkg/dem/SConscript 2009-07-13 15:30:36 UTC (rev 1861)
+++ trunk/pkg/dem/SConscript 2009-07-13 19:31:44 UTC (rev 1862)
@@ -955,7 +955,8 @@
'yade-base',
'yade-multimethods',
'RigidBodyParameters',
- 'Clump'
+ 'Clump',
+ 'VelocityBins'
]),
env.SharedLibrary('SimpleViscoelasticBodyParameters'
Modified: trunk/scripts/test/collider-stride-triax.py
===================================================================
--- trunk/scripts/test/collider-stride-triax.py 2009-07-13 15:30:36 UTC (rev
1861)
+++ trunk/scripts/test/collider-stride-triax.py 2009-07-13 19:31:44 UTC (rev
1862)
@@ -4,9 +4,7 @@
import os.path
loadFrom='/tmp/triax.xml'
if not os.path.exists(loadFrom):
- TriaxialTest(numberOfGrains=8000,fast=True,noFiles=True).load()
- #O.run(500)
- O.save(loadFrom)
+
TriaxialTest(numberOfGrains=20000,fast=True,noFiles=True).generate(loadFrom)
O.load(loadFrom)
log.setLevel('TriaxialCompressionEngine',log.WARN) # shut up
log.setLevel('InsertionSortCollider',log.DEBUG)
@@ -15,9 +13,14 @@
newton=utils.typedEngine('NewtonsDampedLaw')
# use striding; say "if 0:" to disable striding and compare to regular runs
-if 0:
+if 1:
# length by which bboxes will be made larger
collider['sweepLength']=0.05*O.bodies[100].shape['radius']
+ # if this is enabled, bboxes will be enlarged based on velocity bin for
each body
+ if 0:
+ collider['nBins']=3
+ collider['binCoeff']=10
+ log.setLevel('VelocityBins',log.DEBUG)
O.step() # filter out initialization
O.timingEnabled=True
@@ -25,20 +28,12 @@
# run a few times 500 steps, show timings to see what is the trend
# notably, the percentage of collider time should decrease as the max velocity
decreases as well
for i in range(0,5):
- O.run(500,True)
+ O.run(1000,True)
timing.stats()
totalTime+=sum([e.execTime for e in O.engines])
print 'Number of interactions: %d (real ratio:
%g)'%(len(O.interactions),float(O.interactions.countReal())/len(O.interactions))
print '======================================================='
timing.reset()
- # plot velocity histogram
- if 0:
- import pylab
- pylab.hist([sqrt(sum([v**2 for v in b.phys['velocity']])) for b
in O.bodies])
- pylab.title('step %d; max speed
%g'%(O.iter,sqrt(newton['maxVelocitySq'])))
- pylab.grid()
- pylab.savefig('%s-%05d.png'%(loadFrom,O.iter),dpi=40)
- pylab.cla()
print 'Total time: %g s'%(totalTime/1e9)
#quit()
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
