yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07521
[Merge] lp:~bruno-chareyre/yade/collide into lp:yade
Chareyre has proposed merging lp:~bruno-chareyre/yade/collide into lp:yade.
Requested reviews:
Yade developers (yade-dev)
For more details, see:
https://code.launchpad.net/~bruno-chareyre/yade/collide/+merge/59029
Improved colliding (only implemented for non-periodic BCs at the moment).
- Some useless operations are removed (see diff).
- The behaviour with nBins=0 tracks displacements (while nBins>0 tracks max velocity) to update bounds
- Experimental features: targetInterval (adjust bound sizes so that they will be updated each N iterations), updatingDispFactor (only update bounds where the body's motion is at least boundSize/factor, hence saving sorting time).
Overall CPU time reduced by a factor 0.6 (average).
--
https://code.launchpad.net/~bruno-chareyre/yade/collide/+merge/59029
Your team Yade developers is requested to review the proposed merge of lp:~bruno-chareyre/yade/collide into lp:yade.
=== modified file 'core/Bound.hpp'
--- core/Bound.hpp 2010-11-30 13:51:41 +0000
+++ core/Bound.hpp 2011-04-26 08:47:23 +0000
@@ -22,6 +22,13 @@
class Bound: public Serializable, public Indexable{
public:
YADE_CLASS_BASE_DOC_ATTRS_DEPREC_INIT_CTOR_PY(Bound,Serializable,"Object bounding part of space taken by associated body; might be larger, used to optimalize collision detection",
+ #define ORI_VERLET
+ #ifdef ORI_VERLET
+ ((int,lastUpdateIter,0,Attr::readonly,"record iteration of last reference position update |yupdate|"))
+ ((Vector3r,refPos,Vector3r(NaN,NaN,NaN),Attr::readonly,"Reference position, updated at current body position each time the bound dispatcher update bounds - only used if Verlet striding is active |yupdate|"))
+ ((bool,isBounding,false,Attr::readonly,"A flag used to tell when the body moves out of bounds - only used if Verlet striding is active |yupdate|"))
+ ((Real,sweepLength,0, Attr::readonly,"The length used to increase the bounding boxe size - only used if Verlet striding is active |yupdate|"))
+ #endif
((Vector3r,color,Vector3r(1,1,1),,"Color for rendering this object"))
((Vector3r,min,Vector3r(NaN,NaN,NaN),(Attr::noSave | Attr::readonly),"Lower corner of box containing this bound (and the :yref:`Body` as well)"))
((Vector3r,max,Vector3r(NaN,NaN,NaN),(Attr::noSave | Attr::readonly),"Lower corner of box containing this bound (and the :yref:`Body` as well)"))
=== modified file 'core/InteractionContainer.hpp'
--- core/InteractionContainer.hpp 2011-02-27 13:54:43 +0000
+++ core/InteractionContainer.hpp 2011-04-26 08:47:23 +0000
@@ -82,6 +82,9 @@
bool insert(const shared_ptr<Interaction>& i);
bool erase(Body::id_t id1,Body::id_t id2);
const shared_ptr<Interaction>& find(Body::id_t id1,Body::id_t id2);
+// bool found(Body::id_t id1,Body::id_t id2);
+ inline bool found(const Body::id_t& id1,const Body::id_t& id2){
+ assert(bodies); return (id1>id2)?(*bodies)[id2]->intrs.count(id1):(*bodies)[id1]->intrs.count(id2);}
// index access
shared_ptr<Interaction>& operator[](size_t id){return linIntrs[id];}
const shared_ptr<Interaction>& operator[](size_t id) const { return linIntrs[id];}
=== modified file 'pkg/common/Dispatching.cpp'
--- pkg/common/Dispatching.cpp 2011-02-20 10:28:40 +0000
+++ pkg/common/Dispatching.cpp 2011-04-26 08:47:23 +0000
@@ -21,12 +21,33 @@
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
+// #pragma omp parallel for
for(int id=0; id<numBodies; id++){
if(!bodies->exists(id)) continue; // don't delete this check - Janek
const shared_ptr<Body>& b=(*bodies)[id];
shared_ptr<Shape>& shape=b->shape;
if(!shape || !b->isBounded()) continue;
+
+ #ifdef ORI_VERLET
+ if(oriVerlet && b->bound) {
+ Real& sweepLength = b->bound->sweepLength;
+ Vector3r disp = b->state->pos-b->bound->refPos;
+ Real dist = max(abs(disp[0]),max(abs(disp[1]),abs(disp[2])));
+ if (updatingDispFactor<=0 || (dist*updatingDispFactor)>sweepLength) {
+ if (targetInterv>=0) {
+ assert(b->bound->lastUpdateIter<scene->iter);
+ Real sweepMult = dist*targetInterv/((scene->iter-b->bound->lastUpdateIter)*sweepLength);
+ sweepMult = max(0.9,sweepMult);
+ sweepLength=max(0.1*sweepDist,min(sweepLength*sweepMult,sweepDist));
+ } else {
+ sweepLength=sweepDist;
+ }
+ b->bound->isBounding=false;
+ }
+ else continue;
+ }
+ #endif
+
#ifdef BV_FUNCTOR_CACHE
if(!shape->boundFunctor){ shape->boundFunctor=this->getFunctor1D(shape); if(!shape->boundFunctor) continue; }
// LOG_DEBUG("shape->boundFunctor.get()=="<<shape->boundFunctor.get()<<" for "<<b->shape->getClassName()<<", #"<<id);
@@ -37,9 +58,18 @@
#endif
if(!b->bound) continue; // the functor did not create new bound
if(sweepDist>0){
+ #ifdef ORI_VERLET
+ b->bound->refPos=b->state->pos;
+ b->bound->lastUpdateIter=scene->iter;
+ const Real& sweepLength = b->bound->sweepLength;
+ Aabb* aabb=YADE_CAST<Aabb*>(b->bound.get());
+ aabb->min-=Vector3r(sweepLength,sweepLength,sweepLength);
+ aabb->max+=Vector3r(sweepLength,sweepLength,sweepLength);
+ #else
Aabb* aabb=YADE_CAST<Aabb*>(b->bound.get());
aabb->min-=Vector3r(sweepDist,sweepDist,sweepDist);
aabb->max+=Vector3r(sweepDist,sweepDist,sweepDist);
+ #endif
}
if(haveBins){
Aabb* aabb=YADE_CAST<Aabb*>(b->bound.get());
=== modified file 'pkg/common/Dispatching.hpp'
--- pkg/common/Dispatching.hpp 2011-01-29 22:47:18 +0000
+++ pkg/common/Dispatching.hpp 2011-04-26 08:47:23 +0000
@@ -89,6 +89,11 @@
/*additional attrs*/
((bool,activated,true,,"Whether the engine is activated (only should be changed by the collider)"))
((Real,sweepDist,0,,"Distance by which enlarge all bounding boxes, to prevent collider from being run at every step (only should be changed by the collider)."))
+ #ifdef ORI_VERLET
+ ((Real,updatingDispFactor,-1,,"see :yref:`InsertionSortCollider::updatingDispFactor` |yupdate|"))
+ ((Real,targetInterv,-1,,"see :yref:`InsertionSortCollider::targetInterv` |yupdate|"))
+ ((bool,oriVerlet,false,,"Compare Verlet distance with displacement instead of velocity (only should be changed by the collider)"))
+ #endif
,/*ctor*/,/*py*/
);
};
=== modified file 'pkg/common/InsertionSortCollider.cpp'
--- pkg/common/InsertionSortCollider.cpp 2011-02-27 13:54:43 +0000
+++ pkg/common/InsertionSortCollider.cpp 2011-04-26 08:47:23 +0000
@@ -18,44 +18,59 @@
YADE_PLUGIN((InsertionSortCollider))
CREATE_LOGGER(InsertionSortCollider);
+// inlined
// return true if bodies bb overlap in all 3 dimensions
-bool InsertionSortCollider::spatialOverlap(Body::id_t id1, Body::id_t id2) const {
- assert(!periodic);
- return
- (minima[3*id1+0]<=maxima[3*id2+0]) && (maxima[3*id1+0]>=minima[3*id2+0]) &&
- (minima[3*id1+1]<=maxima[3*id2+1]) && (maxima[3*id1+1]>=minima[3*id2+1]) &&
- (minima[3*id1+2]<=maxima[3*id2+2]) && (maxima[3*id1+2]>=minima[3*id2+2]);
-}
+// bool InsertionSortCollider::spatialOverlap(Body::id_t id1, Body::id_t id2) const {
+// assert(!periodic);
+// return
+// (minima[3*id1+0]<=maxima[3*id2+0]) && (maxima[3*id1+0]>=minima[3*id2+0]) &&
+// (minima[3*id1+1]<=maxima[3*id2+1]) && (maxima[3*id1+1]>=minima[3*id2+1]) &&
+// (minima[3*id1+2]<=maxima[3*id2+2]) && (maxima[3*id1+2]>=minima[3*id2+2]);
+// }
// called by the insertion sort if 2 bodies swapped their bounds
void InsertionSortCollider::handleBoundInversion(Body::id_t id1, Body::id_t id2, InteractionContainer* interactions, Scene*){
assert(!periodic);
assert(id1!=id2);
// do bboxes overlap in all 3 dimensions?
- bool overlap=spatialOverlap(id1,id2);
- // existing interaction?
- const shared_ptr<Interaction>& I=interactions->find(id1,id2);
- bool hasInter=(bool)I;
- // interaction doesn't exist and shouldn't, or it exists and should
- if(likely(!overlap && !hasInter)) return;
- if(overlap && hasInter){ return; }
- // create interaction if not yet existing
- if(overlap && !hasInter){ // second condition only for readability
- if(!Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get())) return;
- // LOG_TRACE("Creating new interaction #"<<id1<<"+#"<<id2);
- shared_ptr<Interaction> newI=shared_ptr<Interaction>(new Interaction(id1,id2));
- interactions->insert(newI);
- return;
- }
- if(!overlap && hasInter){ if(!I->isReal()) interactions->erase(id1,id2); return; }
- assert(false); // unreachable
+// bool overlap=spatialOverlap(id1,id2);
+// // existing interaction?
+// bool hasInter=interactions->found(id1,id2);
+// // interaction doesn't exist and shouldn't, or it exists and should
+// if(likely(overlap == hasInter)) return;
+// // create interaction if not yet existing
+// if(overlap && !hasInter){ // second condition only for readability
+// if(!Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get())) return;
+// // LOG_TRACE("Creating new interaction #"<<id1<<"+#"<<id2);
+// shared_ptr<Interaction> newI=shared_ptr<Interaction>(new Interaction(id1,id2));
+// interactions->insert(newI);
+// return;
+// }
+// if(!overlap && hasInter){ if(!interactions->find(id1,id2)->isReal()) interactions->erase(id1,id2); return; }
+// assert(false); // unreachable
+
+ //existing interaction?
+ bool hasInter=interactions->found(id1,id2);
+ if (hasInter && !interactions->find(id1,id2)->isReal() && !spatialOverlap(id1,id2)){ interactions->erase(id1,id2); return;}
+ if (!hasInter && spatialOverlap(id1,id2) && Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get()))
+ interactions->insert(shared_ptr<Interaction>(new Interaction(id1,id2)));
}
void InsertionSortCollider::insertionSort(VecBounds& v, InteractionContainer* interactions, Scene*, bool doCollide){
assert(!periodic);
assert(v.size==(long)v.vec.size());
+ #ifdef ORI_VERLET
+ bool foundUpdatedBound=true;
+ #endif
for(long i=0; i<v.size; i++){
- const Bounds viInit=v[i]; long j=i-1; /* cache hasBB; otherwise 1% overall performance hit */ const bool viInitBB=viInit.flags.hasBB;
+ const Bounds viInit=v[i]; long j=i-1; /* cache hasBB; otherwise 1% overall performance hit */ const bool viInitBB=viInit.flags.hasBB; const bool isMin=viInit.flags.isMin;
+ #ifdef ORI_VERLET
+ if (verletDist>0 && nBins<=0 && viInitBB){
+ if (viInit.isBounding && !foundUpdatedBound) continue;
+ else foundUpdatedBound=true;
+ }
+ #endif
+
while(j>=0 && v[j]>viInit){
v[j+1]=v[j];
#ifdef PISC_DEBUG
@@ -65,10 +80,16 @@
// no collisions without bounding boxes
// also, do not collide body with itself; it sometimes happens for facets aligned perpendicular to an axis, for reasons that are not very clear
// see https://bugs.launchpad.net/yade/+bug/669095
- if(likely(doCollide && viInitBB && v[j].flags.hasBB && (viInit.id!=v[j].id))) handleBoundInversion(viInit.id,v[j].id,interactions,scene);
+ // skip bounds with same isMin flags, since inversion doesn't imply anything in that case
+ if(likely(doCollide && viInitBB && v[j].flags.hasBB && isMin!=v[j].flags.isMin && (viInit.id!=v[j].id))) handleBoundInversion(viInit.id,v[j].id,interactions,scene);
j--;
}
+ #ifdef ORI_VERLET
+ if (j==(i-1) && viInit.isBounding) foundUpdatedBound=false;
+ #endif
v[j+1]=viInit;
+
+
}
}
@@ -94,6 +115,7 @@
}
// STRIDE
+
bool InsertionSortCollider::isActivated(){
// 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
@@ -102,8 +124,18 @@
if(nBins>=1 && newton->velocityBins->checkSize_incrementDists_shouldCollide(scene)) return true;
if(nBins<=0){
if(fastestBodyMaxDist<0){fastestBodyMaxDist=0; return true;}
+ #ifdef ORI_VERLET
+ if (!oriVerlet){
+ fastestBodyMaxDist+=sqrt(newton->maxVelocitySq)*scene->dt;
+ if(fastestBodyMaxDist>=verletDist) return true;
+ } else {
+ fastestBodyMaxDist=newton->maxVelocitySq;
+ if(fastestBodyMaxDist>=1 || fastestBodyMaxDist==0) return true;}
+ #else
fastestBodyMaxDist+=sqrt(newton->maxVelocitySq)*scene->dt;
if(fastestBodyMaxDist>=verletDist) return true;
+ #endif
+
}
if((size_t)BB[0].size!=2*scene->bodies->size()) return true;
if(scene->interactions->dirty) return true;
@@ -159,16 +191,6 @@
// compatibility block, can be removed later
findBoundDispatcherInEnginesIfNoFunctorsAndWarn();
- // update bounds via boundDispatcher
- boundDispatcher->scene=scene;
- boundDispatcher->action();
-
- // if interactions are dirty, force reinitialization
- if(scene->interactions->dirty){
- doInitSort=true;
- scene->interactions->dirty=false;
- }
-
if(verletDist<0){
Real minR=std::numeric_limits<Real>::infinity();
FOREACH(const shared_ptr<Body>& b, *scene->bodies){
@@ -181,15 +203,32 @@
verletDist=isinf(minR) ? 0 : abs(verletDist)*minR;
}
+ // update bounds via boundDispatcher
+ boundDispatcher->scene=scene;
+ #ifdef ORI_VERLET
+ if (oriVerlet) {
+ nBins=0;
+ boundDispatcher->oriVerlet=true;
+ boundDispatcher->sweepDist=verletDist;
+ boundDispatcher->targetInterv=targetInterv;
+ boundDispatcher->updatingDispFactor=updatingDispFactor;}
+ #endif
+ boundDispatcher->action();
+ // if interactions are dirty, force reinitialization
+ if(scene->interactions->dirty){
+ doInitSort=true;
+ scene->interactions->dirty=false;
+ }
+
// STRIDE
- if(verletDist>0){
- // get NewtonIntegrator, to ask for the maximum velocity value
- if(!newton){
- FOREACH(shared_ptr<Engine>& e, scene->engines){ newton=dynamic_pointer_cast<NewtonIntegrator>(e); if(newton) break; }
- if(!newton){ throw runtime_error("InsertionSortCollider.verletDist>0, but unable to locate NewtonIntegrator within O.engines."); }
- }
+ if(verletDist>0){
+ // get NewtonIntegrator, to ask for the maximum velocity value
+ if(!newton){
+ FOREACH(shared_ptr<Engine>& e, scene->engines){ newton=dynamic_pointer_cast<NewtonIntegrator>(e); if(newton) break; }
+ if(!newton){ throw runtime_error("InsertionSortCollider.verletDist>0, but unable to locate NewtonIntegrator within O.engines."); }
}
+ }
ISC_CHECKPOINT("init");
// STRIDE
@@ -203,6 +242,18 @@
// 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.);
+ #ifdef ORI_VERLET
+ newton->oriVerlet=oriVerlet;
+ if (!oriVerlet){
+ Real sweepVelocity=sqrt(newton->maxVelocitySq)*sweepFactor; int stride=-1;
+ if(sweepVelocity>0) {
+ stride=max(1,int((verletDist/sweepVelocity)/scene->dt));
+ boundDispatcher->sweepDist=scene->dt*(stride-1)*sweepVelocity;
+ } else { // no motion
+ boundDispatcher->sweepDist=0; // nothing moves, no need to make bboxes larger
+ }
+ }
+ #else
Real sweepVelocity=sqrt(newton->maxVelocitySq)*sweepFactor; int stride=-1;
if(sweepVelocity>0) {
stride=max(1,int((verletDist/sweepVelocity)/scene->dt));
@@ -212,6 +263,9 @@
}
LOG_DEBUG(scene->time<<"s: stride ≈"<<stride<<"; maxVelocity="<<sqrt(newton->maxVelocitySq)<<", sweepDist="<<boundDispatcher->sweepDist);
fastestBodyMaxDist=0; // reset
+ #endif
+
+
} else { // nBins>=1
if(!newton->velocityBins){ newton->velocityBins=shared_ptr<VelocityBins>(new VelocityBins(nBins,newton->maxVelocitySq,binCoeff,binOverlap)); }
if(!boundDispatcher->velocityBins) boundDispatcher->velocityBins=newton->velocityBins;
@@ -235,9 +289,16 @@
const shared_ptr<Body>& b=Body::byId(id,scene);
if(likely(b)){
const shared_ptr<Bound>& bv=b->bound;
+ #ifdef ORI_VERLET
+ if (oriVerlet){
+ if (bv && bv->isBounding) {BBj[i].isBounding=true; continue;}
+ else BBj[i].isBounding=false;}
+ #endif
// coordinate is min/max if has bounding volume, otherwise both are the position. Add periodic shift so that we are inside the cell
// watch out for the parentheses around ?: within ?: (there was unwanted conversion of the Reals to bools!)
+
BBj[i].coord=((BBj[i].flags.hasBB=((bool)bv)) ? (BBj[i].flags.isMin ? bv->min[j] : bv->max[j]) : (b->state->pos[j])) - (periodic ? BBj.cellDim*BBj[i].period : 0.);
+
} else { BBj[i].flags.hasBB=false; /* for vanished body, keep the coordinate as-is, to minimize inversions. */ }
// if initializing periodic, shift coords & record the period into BBj[i].period
if(doInitSort && periodic) {
@@ -251,6 +312,9 @@
const shared_ptr<Body>& b=Body::byId(id,scene);
if(likely(b)){
const shared_ptr<Bound>& bv=b->bound;
+ #ifdef ORI_VERLET
+ if (oriVerlet && bv && bv->isBounding) continue;
+ #endif
if(likely(bv)) { memcpy(&minima[3*id],&bv->min,3*sizeof(Real)); memcpy(&maxima[3*id],&bv->max,3*sizeof(Real)); } // ⇠faster than 6 assignments
else{ const Vector3r& pos=b->state->pos; memcpy(&minima[3*id],&pos,3*sizeof(Real)); memcpy(&maxima[3*id],&pos,3*sizeof(Real)); }
} else { memset(&minima[3*id],0,3*sizeof(Real)); memset(&maxima[3*id],0,3*sizeof(Real)); }
=== modified file 'pkg/common/InsertionSortCollider.hpp'
--- pkg/common/InsertionSortCollider.hpp 2011-01-09 16:34:50 +0000
+++ pkg/common/InsertionSortCollider.hpp 2011-04-26 08:47:23 +0000
@@ -68,7 +68,7 @@
// #define this macro to enable timing within this engine
-//#define ISC_TIMING
+// #define ISC_TIMING
// #define to turn on some tracing information for the periodic part
// all code under this can be probably removed at some point, when the collider will have been tested thoroughly
@@ -90,6 +90,9 @@
Real coord;
//! id of the body this bound belongs to
Body::id_t id;
+ #ifdef ORI_VERLET
+ bool isBounding;
+ #endif
//! periodic cell coordinate
int period;
//! is it the minimum (true) or maximum (false) bound?
@@ -152,12 +155,19 @@
*/
void insertionSort(VecBounds& v,InteractionContainer*,Scene*,bool doCollide=true);
void handleBoundInversion(Body::id_t,Body::id_t,InteractionContainer*,Scene*);
- bool spatialOverlap(Body::id_t,Body::id_t) const;
+// bool spatialOverlap(Body::id_t,Body::id_t) const;
// periodic variants
void insertionSortPeri(VecBounds& v,InteractionContainer*,Scene*,bool doCollide=true);
void handleBoundInversionPeri(Body::id_t,Body::id_t,InteractionContainer*,Scene*);
bool spatialOverlapPeri(Body::id_t,Body::id_t,Scene*,Vector3i&) const;
+ inline bool spatialOverlap(const Body::id_t& id1, const Body::id_t& id2) const {
+ assert(!periodic);
+ return (minima[3*id1+0]<=maxima[3*id2+0]) && (maxima[3*id1+0]>=minima[3*id2+0]) &&
+ (minima[3*id1+1]<=maxima[3*id2+1]) && (maxima[3*id1+1]>=minima[3*id2+1]) &&
+ (minima[3*id1+2]<=maxima[3*id2+2]) && (maxima[3*id1+2]>=minima[3*id2+2]);
+}
+
static Real cellWrap(const Real, const Real, const Real, int&);
static Real cellWrapRel(const Real, const Real, const Real);
@@ -199,7 +209,12 @@
",
((int,sortAxis,0,,"Axis for the initial contact detection."))
((bool,sortThenCollide,false,,"Separate sorting and colliding phase; it is MUCH slower, but all interactions are processed at every step; this effectively makes the collider non-persistent, not remembering last state. (The default behavior relies on the fact that inversions during insertion sort are overlaps of bounding boxes that just started/ceased to exist, and only processes those; this makes the collider much more efficient.)"))
- ((Real,verletDist,((void)"Automatically initialized",-.05),,"Length by which to enlarge particle bounds, to avoid running collider at every step. Stride disabled if zero. Negative value will trigger automatic computation, so that the real value will be |verletDist| × minimum spherical particle radius; if there are no spherical particles, it will be disabled."))
+ #ifdef ORI_VERLET
+ ((bool,oriVerlet,true,,"Compare Verlet distance with displacement instead of velocity (only used if nBins<=0)"))
+ ((int,targetInterv,-1,,"(experimental) Target number of iterations between bound update, used to define a smaller sweep distance for slower grains if >0, else always use 1*verletDist. Only used if oriVerlet=true."))
+ ((Real,updatingDispFactor,-1,,"(experimental) Displacement factor used to trigger bound update when oriVerlet=true: the bound is updated only if updatingDispFactor*disp>sweepDist when >0, else all bounds are updated."))
+ #endif
+ ((Real,verletDist,((void)"Automatically initialized",-.15),,"Length by which to enlarge particle bounds, to avoid running collider at every step. Stride disabled if zero. Negative value will trigger automatic computation, so that the real value will be |verletDist| × minimum spherical particle radius; if there are no spherical particles, it will be disabled."))
((Real,sweepFactor,1.05,,"Overestimation factor for the sweep velocity; must be >=1.0. Has no influence on verletDist, only on the computed stride. [DEPRECATED, is used only when bins are not used]."))
((Real,fastestBodyMaxDist,-1,,"Maximum displacement of the fastest body since last run; if >= verletDist, we could get out of bboxes and will trigger full run. DEPRECATED, was only used without bins. |yupdate|"))
((int,nBins,5,,"Number of velocity bins for striding. If <=0, bin-less strigin is used (this is however DEPRECATED)."))
=== modified file 'pkg/common/InteractionLoop.cpp'
--- pkg/common/InteractionLoop.cpp 2011-02-27 14:26:15 +0000
+++ pkg/common/InteractionLoop.cpp 2011-04-26 08:47:23 +0000
@@ -103,7 +103,11 @@
bool swap=false;
// IGeomDispatcher
+ #ifdef ORI_VERLET
+ if(unlikely(!I->functorCache.geom)){
+ #else
if(unlikely(!I->functorCache.geom || !I->functorCache.phys)){
+ #endif
I->functorCache.geom=geomDispatcher->getFunctor2D(b1_->shape,b2_->shape,swap);
// returns NULL ptr if no functor exists; remember that and shortcut
if(!I->functorCache.geom) {I->functorCache.geomExists=false; continue; }
@@ -112,10 +116,10 @@
// arguments for the geom functor are in the reverse order (dispatcher would normally call goReverse).
// we don't remember the fact that is reverse, so we swap bodies within the interaction
// and can call go in all cases
- if(unlikely(swap)){I->swapOrder(); }
+ if(unlikely(swap)){I->swapOrder();}
// body pointers must be updated, in case we swapped
- const shared_ptr<Body>& b1=Body::byId(I->getId1(),scene);
- const shared_ptr<Body>& b2=Body::byId(I->getId2(),scene);
+ const shared_ptr<Body>& b1=swap?b2_:b1_;
+ const shared_ptr<Body>& b2=swap?b1_:b2_;
assert(I->functorCache.geom);
bool wasReal=I->isReal();
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2011-04-22 09:01:59 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2011-04-26 08:47:23 +0000
@@ -9,7 +9,6 @@
#include<yade/pkg/dem/NewtonIntegrator.hpp>
#include<yade/core/Scene.hpp>
#include<yade/pkg/dem/Clump.hpp>
-#include<yade/pkg/common/VelocityBins.hpp>
#include<yade/lib/base/Math.hpp>
YADE_PLUGIN((NewtonIntegrator));
@@ -66,7 +65,19 @@
maxVelocitySq=max(maxVelocitySq,state->vel.squaredNorm());
#endif
}
-
+#ifdef ORI_VERLET
+void NewtonIntegrator::saveMaximaDisplacement(const shared_ptr<Body>& b){
+ Vector3r disp=b->state->pos-b->bound->refPos;
+ Real maxDisp=max(abs(disp[0]),max(abs(disp[1]),abs(disp[2])));//b->bound->sweepLength;
+ if (maxDisp<=b->bound->sweepLength) b->bound->isBounding = true; else b->bound->isBounding = false;
+ if (b->bound->sweepLength>0) maxDisp/=b->bound->sweepLength; else maxDisp=(maxDisp>0)?2:0;
+ #ifdef YADE_OPENMP
+ Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,maxDisp);
+ #else
+ maxVelocitySq=max(maxVelocitySq,maxDisp);
+ #endif
+}
+#endif
void NewtonIntegrator::action()
{
scene->forces.sync();
@@ -108,21 +119,7 @@
if(unlikely(!b || b->isClumpMember())) continue;
State* state=b->state.get(); const Body::id_t& id=b->getId();
- Vector3r f=Vector3r::Zero(), m=Vector3r::Zero();
- // clumps forces
- if(b->isClump()) {
- b->shape->cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
- #ifdef YADE_OPENMP
- //it is safe here, since onky one thread will read/write
- scene->forces.getTorqueUnsynced(id)+=m;
- scene->forces.getForceUnsynced(id)+=f;
- #else
- scene->forces.addTorque(id,m);
- scene->forces.addForce(id, f);
- #endif
- }
- //in most cases, the initial force on clumps will zero and next line is not changing f and m, but make sure we don't miss something (e.g. user defined forces on clumps)
- f=scene->forces.getForce(id); m=scene->forces.getTorque(id);
+ Vector3r f=scene->forces.getForce(id), m=scene->forces.getTorque(id);
#ifdef YADE_DEBUG
if(isnan(f[0])||isnan(f[1])||isnan(f[2])) throw runtime_error(("NewtonIntegrator: NaN force acting on #"+lexical_cast<string>(id)+".").c_str());
if(isnan(m[0])||isnan(m[1])||isnan(m[2])) throw runtime_error(("NewtonIntegrator: NaN torque acting on #"+lexical_cast<string>(id)+".").c_str());
@@ -134,8 +131,6 @@
// in aperiodic boundaries, it is equal to absolute velocity
Vector3r fluctVel=isPeriodic?scene->cell->bodyFluctuationVel(b->state->pos,b->state->vel):state->vel;
-
-
// numerical damping & kinetic energy
if(unlikely(trackEnergy)) updateEnergy(b,state,fluctVel,f,m);
@@ -144,6 +139,8 @@
// for particles not totally blocked, compute accelerations; otherwise, the computations would be useless
if (state->blockedDOFs!=State::DOF_ALL) {
+ // forces
+ if(b->isClump()) b->shape->cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
// linear acceleration
Vector3r linAccel=computeAccel(f,state->mass,state->blockedDOFs);
cundallDamp2nd(dt,f,fluctVel,linAccel);
@@ -163,7 +160,10 @@
leapfrogTranslate(state,id,dt);
if(!useAspherical) leapfrogSphericalRotate(state,id,dt);
else leapfrogAsphericalRotate(state,id,dt,m);
-
+ #ifdef ORI_VERLET
+ if (oriVerlet) saveMaximaDisplacement(b);
+ else
+ #endif
saveMaximaVelocity(id,state);
// move individual members of the clump, save maxima velocity (for collider stride)
if(b->isClump()) Clump::moveMembers(b,scene,this);
=== modified file 'pkg/dem/NewtonIntegrator.hpp'
--- pkg/dem/NewtonIntegrator.hpp 2011-04-22 09:01:59 +0000
+++ pkg/dem/NewtonIntegrator.hpp 2011-04-26 08:47:23 +0000
@@ -11,6 +11,7 @@
#include<yade/core/Interaction.hpp>
#include<yade/lib/base/Math.hpp>
#include<yade/pkg/common/Callbacks.hpp>
+#include<yade/pkg/common/VelocityBins.hpp>
#ifdef YADE_OPENMP
#include<omp.h>
#endif
@@ -46,6 +47,9 @@
public:
// function to save maximum velocity, for the verlet-distance optimization
void saveMaximaVelocity(const Body::id_t& id, State* state);
+ #ifdef ORI_VERLET
+ void saveMaximaDisplacement(const shared_ptr<Body>& b);
+ #endif
#ifdef YADE_OPENMP
vector<Real> threadMaxVelocitySq;
#endif
@@ -56,6 +60,7 @@
((Real,damping,0.2,,"damping coefficient for Cundall's non viscous damping (see [Chareyre2005]_) [-]"))
((Real,maxVelocitySq,NaN,,"store square of max. velocity, for informative purposes; computed again at every step. |yupdate|"))
((bool,exactAsphericalRot,true,,"Enable more exact body rotation integrator for :yref:`aspherical bodies<Body.aspherical>` *only*, using formulation from [Allen1989]_, pg. 89."))
+ ((bool,oriVerlet,false,,"Compare Verlet distance with displacement instead of velocity (updated by the collider)"))
((Matrix3r,prevVelGrad,Matrix3r::Zero(),,"Store previous velocity gradient (:yref:`Cell::velGrad`) to track acceleration. |yupdate|"))
#ifdef YADE_BODY_CALLBACK
((vector<shared_ptr<BodyCallback> >,callbacks,,,"List (std::vector in c++) of :yref:`BodyCallbacks<BodyCallback>` which will be called for each body as it is being processed."))
Follow ups
