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