yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #07522
Re: [Merge] lp:~bruno-chareyre/yade/collide into lp:yade
Hi, Bruno!
Please, do not hurry up with merging. Let's test it.
Thanks.
Anton
On Tue, Apr 26, 2011 at 10:47 AM, Chareyre
<bruno.chareyre@xxxxxxxxxxxxxxx> wrote:
> 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."))
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>
>
Follow ups
References
