Yade

Merge lp:~bruno-chareyre/yade/collide into lp:~yade-dev/yade/trunk

  1. collide
  2. Merge into trunk
Proposed by Bruno Chareyre
Status: Merged
Merged at revision: 2922
Proposed branch: lp:~bruno-chareyre/yade/collide
Merge into: lp:~yade-dev/yade/trunk
Diff against target: 552 lines (+194/-71)
9 files modified
core/Bound.hpp (+7/-0)
core/InteractionContainer.hpp (+3/-0)
pkg/common/Dispatching.cpp (+32/-2)
pkg/common/Dispatching.hpp (+6/-0)
pkg/common/InsertionSortCollider.cpp (+88/-43)
pkg/common/InsertionSortCollider.hpp (+22/-3)
pkg/common/InteractionLoop.cpp (+7/-3)
pkg/dem/NewtonIntegrator.cpp (+22/-20)
pkg/dem/NewtonIntegrator.hpp (+7/-0)
To merge this branch: bzr merge lp:~bruno-chareyre/yade/collide
Reviewer Review Type Date Requested Status
Yade developers Pending
Review via email: mp+59029@code.launchpad.net

Description of the change

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).

To post a comment you must log in.
lp:~bruno-chareyre/yade/collide updated
2823. By Bruno Chareyre

- default params + a smarter inversion handling for not-overlapping bounds

2824. By Bruno Chareyre

-no change

2825. By Bruno Chareyre

-merge

2826. By Bruno Chareyre

- fix commit mistake + documentation

2827. By Bruno Chareyre

- missing bracket

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'core/Bound.hpp'
2--- core/Bound.hpp 2010-11-30 13:51:41 +0000
3+++ core/Bound.hpp 2011-04-26 10:44:50 +0000
4@@ -22,6 +22,13 @@
5 class Bound: public Serializable, public Indexable{
6 public:
7 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",
8+ #define ORI_VERLET
9+ #ifdef ORI_VERLET
10+ ((int,lastUpdateIter,0,Attr::readonly,"record iteration of last reference position update |yupdate|"))
11+ ((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 oriVerlet striding is active |yupdate|"))
12+ ((bool,isBounding,false,Attr::readonly,"A flag used to tell when the body moves out of bounds - only used if oriVerlet striding is active :yref:`BoundDispatcher::updatingDispFactor`>0 |yupdate|"))
13+ ((Real,sweepLength,0, Attr::readonly,"The length used to increase the bounding boxe size, can be adjusted on the basis of previous displacement if :yref:`BoundDispatcher::targetInterv`>0. Only used if oriVerlet striding is active (:yref:`InsertionSortCollider::oriVerlet`=true) |yupdate|"))
14+ #endif
15 ((Vector3r,color,Vector3r(1,1,1),,"Color for rendering this object"))
16 ((Vector3r,min,Vector3r(NaN,NaN,NaN),(Attr::noSave | Attr::readonly),"Lower corner of box containing this bound (and the :yref:`Body` as well)"))
17 ((Vector3r,max,Vector3r(NaN,NaN,NaN),(Attr::noSave | Attr::readonly),"Lower corner of box containing this bound (and the :yref:`Body` as well)"))
18
19=== modified file 'core/InteractionContainer.hpp'
20--- core/InteractionContainer.hpp 2011-02-27 13:54:43 +0000
21+++ core/InteractionContainer.hpp 2011-04-26 10:44:50 +0000
22@@ -82,6 +82,9 @@
23 bool insert(const shared_ptr<Interaction>& i);
24 bool erase(Body::id_t id1,Body::id_t id2);
25 const shared_ptr<Interaction>& find(Body::id_t id1,Body::id_t id2);
26+// bool found(Body::id_t id1,Body::id_t id2);
27+ inline bool found(const Body::id_t& id1,const Body::id_t& id2){
28+ assert(bodies); return (id1>id2)?(*bodies)[id2]->intrs.count(id1):(*bodies)[id1]->intrs.count(id2);}
29 // index access
30 shared_ptr<Interaction>& operator[](size_t id){return linIntrs[id];}
31 const shared_ptr<Interaction>& operator[](size_t id) const { return linIntrs[id];}
32
33=== modified file 'pkg/common/Dispatching.cpp'
34--- pkg/common/Dispatching.cpp 2011-02-20 10:28:40 +0000
35+++ pkg/common/Dispatching.cpp 2011-04-26 10:44:50 +0000
36@@ -21,12 +21,33 @@
37 const long numBodies=(long)bodies->size();
38 bool haveBins=(bool)velocityBins;
39 if(sweepDist!=0 && haveBins){ LOG_FATAL("Only one of sweepDist or velocityBins can used!"); abort(); }
40- //#pragma omp parallel for
41+// #pragma omp parallel for
42 for(int id=0; id<numBodies; id++){
43 if(!bodies->exists(id)) continue; // don't delete this check - Janek
44 const shared_ptr<Body>& b=(*bodies)[id];
45 shared_ptr<Shape>& shape=b->shape;
46- if(!shape || !b->isBounded()) continue;
47+ if(!shape || !b->isBounded()) continue;
48+ #ifdef ORI_VERLET
49+ if(oriVerlet && b->bound) {
50+ //Do we need to update this bound (determined in Newton)?
51+ if(!b->bound->isBounding){
52+ Real& sweepLength = b->bound->sweepLength;
53+ if (targetInterv>=0) {
54+ Vector3r disp = b->state->pos-b->bound->refPos;
55+ Real dist = max(abs(disp[0]),max(abs(disp[1]),abs(disp[2])));
56+ assert(b->bound->lastUpdateIter<scene->iter);
57+// Real sweepMult = dist*targetInterv/((scene->iter-b->bound->lastUpdateIter)*sweepLength);
58+ Real newLength = dist*targetInterv/(scene->iter-b->bound->lastUpdateIter);
59+ newLength = max(0.9*sweepLength,newLength);//don't decrease size too fast to prevent time consuming oscillations
60+ sweepLength=max(minSweepDistFactor*sweepDist,min(newLength,sweepDist));
61+ } else {
62+ sweepLength=sweepDist;
63+ }
64+ }
65+ else continue;
66+ }
67+ #endif
68+
69 #ifdef BV_FUNCTOR_CACHE
70 if(!shape->boundFunctor){ shape->boundFunctor=this->getFunctor1D(shape); if(!shape->boundFunctor) continue; }
71 // LOG_DEBUG("shape->boundFunctor.get()=="<<shape->boundFunctor.get()<<" for "<<b->shape->getClassName()<<", #"<<id);
72@@ -37,9 +58,18 @@
73 #endif
74 if(!b->bound) continue; // the functor did not create new bound
75 if(sweepDist>0){
76+ #ifdef ORI_VERLET
77+ b->bound->refPos=b->state->pos;
78+ b->bound->lastUpdateIter=scene->iter;
79+ const Real& sweepLength = b->bound->sweepLength;
80+ Aabb* aabb=YADE_CAST<Aabb*>(b->bound.get());
81+ aabb->min-=Vector3r(sweepLength,sweepLength,sweepLength);
82+ aabb->max+=Vector3r(sweepLength,sweepLength,sweepLength);
83+ #else
84 Aabb* aabb=YADE_CAST<Aabb*>(b->bound.get());
85 aabb->min-=Vector3r(sweepDist,sweepDist,sweepDist);
86 aabb->max+=Vector3r(sweepDist,sweepDist,sweepDist);
87+ #endif
88 }
89 if(haveBins){
90 Aabb* aabb=YADE_CAST<Aabb*>(b->bound.get());
91
92=== modified file 'pkg/common/Dispatching.hpp'
93--- pkg/common/Dispatching.hpp 2011-01-29 22:47:18 +0000
94+++ pkg/common/Dispatching.hpp 2011-04-26 10:44:50 +0000
95@@ -89,6 +89,12 @@
96 /*additional attrs*/
97 ((bool,activated,true,,"Whether the engine is activated (only should be changed by the collider)"))
98 ((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)."))
99+ #ifdef ORI_VERLET
100+ ((Real,minSweepDistFactor,0.05,,"Minimal distance by which enlarge all bounding boxes; superseeds computed value of sweepDist when lower that minSweepDist."))
101+ ((Real,updatingDispFactor,-1,,"see :yref:`InsertionSortCollider::updatingDispFactor` |yupdate|"))
102+ ((Real,targetInterv,-1,,"see :yref:`InsertionSortCollider::targetInterv` |yupdate|"))
103+ ((bool,oriVerlet,false,,"Compare Verlet distance with displacement instead of velocity (only should be changed by the collider)"))
104+ #endif
105 ,/*ctor*/,/*py*/
106 );
107 };
108
109=== modified file 'pkg/common/InsertionSortCollider.cpp'
110--- pkg/common/InsertionSortCollider.cpp 2011-02-27 13:54:43 +0000
111+++ pkg/common/InsertionSortCollider.cpp 2011-04-26 10:44:50 +0000
112@@ -18,44 +18,39 @@
113 YADE_PLUGIN((InsertionSortCollider))
114 CREATE_LOGGER(InsertionSortCollider);
115
116-// return true if bodies bb overlap in all 3 dimensions
117-bool InsertionSortCollider::spatialOverlap(Body::id_t id1, Body::id_t id2) const {
118- assert(!periodic);
119- return
120- (minima[3*id1+0]<=maxima[3*id2+0]) && (maxima[3*id1+0]>=minima[3*id2+0]) &&
121- (minima[3*id1+1]<=maxima[3*id2+1]) && (maxima[3*id1+1]>=minima[3*id2+1]) &&
122- (minima[3*id1+2]<=maxima[3*id2+2]) && (maxima[3*id1+2]>=minima[3*id2+2]);
123-}
124
125 // called by the insertion sort if 2 bodies swapped their bounds
126 void InsertionSortCollider::handleBoundInversion(Body::id_t id1, Body::id_t id2, InteractionContainer* interactions, Scene*){
127 assert(!periodic);
128 assert(id1!=id2);
129- // do bboxes overlap in all 3 dimensions?
130- bool overlap=spatialOverlap(id1,id2);
131- // existing interaction?
132+ //existing interaction?
133+ if (interactions->found(id1,id2)) {
134+ //if existing interaction is virtual and bounds don't overlap, remove it
135+ if (!interactions->find(id1,id2)->isReal() && !spatialOverlap(id1,id2)){ interactions->erase(id1,id2); return;}}
136+ //if it doesn't exist and bounds overlap, create a virtual interaction
137+ else if (spatialOverlap(id1,id2) && Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get()))
138+ interactions->insert(shared_ptr<Interaction>(new Interaction(id1,id2)));
139+}
140+
141+// called by the insertion sort if 2 bodies swapped their bounds
142+void InsertionSortCollider::handleBoundSplit(Body::id_t id1, Body::id_t id2, InteractionContainer* interactions, Scene*){
143+ assert(!periodic);
144+ assert(id1!=id2);
145 const shared_ptr<Interaction>& I=interactions->find(id1,id2);
146- bool hasInter=(bool)I;
147- // interaction doesn't exist and shouldn't, or it exists and should
148- if(likely(!overlap && !hasInter)) return;
149- if(overlap && hasInter){ return; }
150- // create interaction if not yet existing
151- if(overlap && !hasInter){ // second condition only for readability
152- if(!Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get())) return;
153- // LOG_TRACE("Creating new interaction #"<<id1<<"+#"<<id2);
154- shared_ptr<Interaction> newI=shared_ptr<Interaction>(new Interaction(id1,id2));
155- interactions->insert(newI);
156- return;
157- }
158- if(!overlap && hasInter){ if(!I->isReal()) interactions->erase(id1,id2); return; }
159- assert(false); // unreachable
160+ if ((bool)I && !I->isReal()) interactions->erase(id1,id2);
161 }
162
163 void InsertionSortCollider::insertionSort(VecBounds& v, InteractionContainer* interactions, Scene*, bool doCollide){
164 assert(!periodic);
165 assert(v.size==(long)v.vec.size());
166- for(long i=0; i<v.size; i++){
167+ for(long i=1; i<v.size; i++){
168+ #ifdef ORI_VERLET
169+ assert (v[i].flags.hasBB && v[i-1].flags.hasBB);
170+ if (updatingDispFactor>0 && v[i].isBounding && v[i-1].isBounding) continue;
171+ #endif
172 const Bounds viInit=v[i]; long j=i-1; /* cache hasBB; otherwise 1% overall performance hit */ const bool viInitBB=viInit.flags.hasBB;
173+ const bool isMin=viInit.flags.isMin;
174+
175 while(j>=0 && v[j]>viInit){
176 v[j+1]=v[j];
177 #ifdef PISC_DEBUG
178@@ -65,7 +60,11 @@
179 // no collisions without bounding boxes
180 // also, do not collide body with itself; it sometimes happens for facets aligned perpendicular to an axis, for reasons that are not very clear
181 // see https://bugs.launchpad.net/yade/+bug/669095
182- if(likely(doCollide && viInitBB && v[j].flags.hasBB && (viInit.id!=v[j].id))) handleBoundInversion(viInit.id,v[j].id,interactions,scene);
183+ // skip bounds with same isMin flags, since inversion doesn't imply anything in that case
184+ if(isMin!=v[j].flags.isMin && likely(doCollide && viInitBB && v[j].flags.hasBB && (viInit.id!=v[j].id))) {
185+ if (isMin) handleBoundInversion(viInit.id,v[j].id,interactions,scene);
186+ else handleBoundSplit(viInit.id,v[j].id,interactions,scene);
187+ }
188 j--;
189 }
190 v[j+1]=viInit;
191@@ -94,6 +93,7 @@
192 }
193
194 // STRIDE
195+
196 bool InsertionSortCollider::isActivated(){
197 // activated if number of bodies changes (hence need to refresh collision information)
198 // or the time of scheduled run already came, or we were never scheduled yet
199@@ -102,8 +102,18 @@
200 if(nBins>=1 && newton->velocityBins->checkSize_incrementDists_shouldCollide(scene)) return true;
201 if(nBins<=0){
202 if(fastestBodyMaxDist<0){fastestBodyMaxDist=0; return true;}
203+ #ifdef ORI_VERLET
204+ if (!oriVerlet){
205+ fastestBodyMaxDist+=sqrt(newton->maxVelocitySq)*scene->dt;
206+ if(fastestBodyMaxDist>=verletDist) return true;
207+ } else {
208+ fastestBodyMaxDist=newton->maxVelocitySq;
209+ if(fastestBodyMaxDist>=1 || fastestBodyMaxDist==0) return true;}
210+ #else
211 fastestBodyMaxDist+=sqrt(newton->maxVelocitySq)*scene->dt;
212 if(fastestBodyMaxDist>=verletDist) return true;
213+ #endif
214+
215 }
216 if((size_t)BB[0].size!=2*scene->bodies->size()) return true;
217 if(scene->interactions->dirty) return true;
218@@ -159,16 +169,6 @@
219 // compatibility block, can be removed later
220 findBoundDispatcherInEnginesIfNoFunctorsAndWarn();
221
222- // update bounds via boundDispatcher
223- boundDispatcher->scene=scene;
224- boundDispatcher->action();
225-
226- // if interactions are dirty, force reinitialization
227- if(scene->interactions->dirty){
228- doInitSort=true;
229- scene->interactions->dirty=false;
230- }
231-
232 if(verletDist<0){
233 Real minR=std::numeric_limits<Real>::infinity();
234 FOREACH(const shared_ptr<Body>& b, *scene->bodies){
235@@ -181,15 +181,33 @@
236 verletDist=isinf(minR) ? 0 : abs(verletDist)*minR;
237 }
238
239+ // update bounds via boundDispatcher
240+ boundDispatcher->scene=scene;
241+ #ifdef ORI_VERLET
242+ if (oriVerlet) {
243+ nBins=0;
244+ if(boundDispatcher->velocityBins) boundDispatcher->velocityBins=shared_ptr<VelocityBins>();//avoid FATAL in boundDispatcher->action() when switching from nBin>0 to oriVerlet=true
245+ boundDispatcher->oriVerlet=true;
246+ boundDispatcher->sweepDist=verletDist;
247+ boundDispatcher->targetInterv=targetInterv;
248+ boundDispatcher->updatingDispFactor=updatingDispFactor;}
249+ #endif
250+ boundDispatcher->action();
251
252+ // if interactions are dirty, force reinitialization
253+ if(scene->interactions->dirty){
254+ doInitSort=true;
255+ scene->interactions->dirty=false;
256+ }
257+
258 // STRIDE
259- if(verletDist>0){
260- // get NewtonIntegrator, to ask for the maximum velocity value
261- if(!newton){
262- FOREACH(shared_ptr<Engine>& e, scene->engines){ newton=dynamic_pointer_cast<NewtonIntegrator>(e); if(newton) break; }
263- if(!newton){ throw runtime_error("InsertionSortCollider.verletDist>0, but unable to locate NewtonIntegrator within O.engines."); }
264- }
265+ if(verletDist>0){
266+ // get NewtonIntegrator, to ask for the maximum velocity value
267+ if(!newton){
268+ FOREACH(shared_ptr<Engine>& e, scene->engines){ newton=dynamic_pointer_cast<NewtonIntegrator>(e); if(newton) break; }
269+ if(!newton){ throw runtime_error("InsertionSortCollider.verletDist>0, but unable to locate NewtonIntegrator within O.engines."); }
270 }
271+ }
272 ISC_CHECKPOINT("init");
273
274 // STRIDE
275@@ -203,6 +221,20 @@
276 // reset bins, in case they were active but are not anymore
277 if(newton->velocityBins) newton->velocityBins=shared_ptr<VelocityBins>(); if(boundDispatcher->velocityBins) boundDispatcher->velocityBins=shared_ptr<VelocityBins>();
278 assert(strideActive); assert(newton->maxVelocitySq>=0); assert(sweepFactor>1.);
279+ #ifdef ORI_VERLET
280+ newton->oriVerlet=oriVerlet;
281+ newton->updatingDispFactor=updatingDispFactor;
282+
283+ if (!oriVerlet){
284+ Real sweepVelocity=sqrt(newton->maxVelocitySq)*sweepFactor; int stride=-1;
285+ if(sweepVelocity>0) {
286+ stride=max(1,int((verletDist/sweepVelocity)/scene->dt));
287+ boundDispatcher->sweepDist=scene->dt*(stride-1)*sweepVelocity;
288+ } else { // no motion
289+ boundDispatcher->sweepDist=0; // nothing moves, no need to make bboxes larger
290+ }
291+ }
292+ #else
293 Real sweepVelocity=sqrt(newton->maxVelocitySq)*sweepFactor; int stride=-1;
294 if(sweepVelocity>0) {
295 stride=max(1,int((verletDist/sweepVelocity)/scene->dt));
296@@ -212,6 +244,9 @@
297 }
298 LOG_DEBUG(scene->time<<"s: stride ≈"<<stride<<"; maxVelocity="<<sqrt(newton->maxVelocitySq)<<", sweepDist="<<boundDispatcher->sweepDist);
299 fastestBodyMaxDist=0; // reset
300+ #endif
301+
302+
303 } else { // nBins>=1
304 if(!newton->velocityBins){ newton->velocityBins=shared_ptr<VelocityBins>(new VelocityBins(nBins,newton->maxVelocitySq,binCoeff,binOverlap)); }
305 if(!boundDispatcher->velocityBins) boundDispatcher->velocityBins=newton->velocityBins;
306@@ -235,9 +270,16 @@
307 const shared_ptr<Body>& b=Body::byId(id,scene);
308 if(likely(b)){
309 const shared_ptr<Bound>& bv=b->bound;
310+ #ifdef ORI_VERLET
311+ if (oriVerlet){
312+ if (bv && bv->isBounding) {BBj[i].isBounding=true; continue;}
313+ else BBj[i].isBounding=false;}
314+ #endif
315 // coordinate is min/max if has bounding volume, otherwise both are the position. Add periodic shift so that we are inside the cell
316 // watch out for the parentheses around ?: within ?: (there was unwanted conversion of the Reals to bools!)
317+
318 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.);
319+
320 } else { BBj[i].flags.hasBB=false; /* for vanished body, keep the coordinate as-is, to minimize inversions. */ }
321 // if initializing periodic, shift coords & record the period into BBj[i].period
322 if(doInitSort && periodic) {
323@@ -251,6 +293,9 @@
324 const shared_ptr<Body>& b=Body::byId(id,scene);
325 if(likely(b)){
326 const shared_ptr<Bound>& bv=b->bound;
327+ #ifdef ORI_VERLET
328+ if (oriVerlet && bv && bv->isBounding) continue;
329+ #endif
330 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
331 else{ const Vector3r& pos=b->state->pos; memcpy(&minima[3*id],&pos,3*sizeof(Real)); memcpy(&maxima[3*id],&pos,3*sizeof(Real)); }
332 } else { memset(&minima[3*id],0,3*sizeof(Real)); memset(&maxima[3*id],0,3*sizeof(Real)); }
333
334=== modified file 'pkg/common/InsertionSortCollider.hpp'
335--- pkg/common/InsertionSortCollider.hpp 2011-01-09 16:34:50 +0000
336+++ pkg/common/InsertionSortCollider.hpp 2011-04-26 10:44:50 +0000
337@@ -68,7 +68,7 @@
338
339
340 // #define this macro to enable timing within this engine
341-//#define ISC_TIMING
342+// #define ISC_TIMING
343
344 // #define to turn on some tracing information for the periodic part
345 // all code under this can be probably removed at some point, when the collider will have been tested thoroughly
346@@ -90,6 +90,9 @@
347 Real coord;
348 //! id of the body this bound belongs to
349 Body::id_t id;
350+ #ifdef ORI_VERLET
351+ bool isBounding;
352+ #endif
353 //! periodic cell coordinate
354 int period;
355 //! is it the minimum (true) or maximum (false) bound?
356@@ -152,12 +155,21 @@
357 */
358 void insertionSort(VecBounds& v,InteractionContainer*,Scene*,bool doCollide=true);
359 void handleBoundInversion(Body::id_t,Body::id_t,InteractionContainer*,Scene*);
360- bool spatialOverlap(Body::id_t,Body::id_t) const;
361+// bool spatialOverlap(Body::id_t,Body::id_t) const;
362
363 // periodic variants
364 void insertionSortPeri(VecBounds& v,InteractionContainer*,Scene*,bool doCollide=true);
365 void handleBoundInversionPeri(Body::id_t,Body::id_t,InteractionContainer*,Scene*);
366+ void handleBoundSplit(Body::id_t,Body::id_t,InteractionContainer*,Scene*);
367+
368 bool spatialOverlapPeri(Body::id_t,Body::id_t,Scene*,Vector3i&) const;
369+ inline bool spatialOverlap(const Body::id_t& id1, const Body::id_t& id2) const {
370+ assert(!periodic);
371+ return (minima[3*id1+0]<=maxima[3*id2+0]) && (maxima[3*id1+0]>=minima[3*id2+0]) &&
372+ (minima[3*id1+1]<=maxima[3*id2+1]) && (maxima[3*id1+1]>=minima[3*id2+1]) &&
373+ (minima[3*id1+2]<=maxima[3*id2+2]) && (maxima[3*id1+2]>=minima[3*id2+2]);
374+ }
375+
376 static Real cellWrap(const Real, const Real, const Real, int&);
377 static Real cellWrapRel(const Real, const Real, const Real);
378
379@@ -199,7 +211,14 @@
380 ",
381 ((int,sortAxis,0,,"Axis for the initial contact detection."))
382 ((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.)"))
383- ((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."))
384+ #ifdef ORI_VERLET
385+ ((bool,oriVerlet,true,,"Compare Verlet distance with displacement instead of velocity (only used if nBins<=0)"))
386+ ((int,targetInterv,30,,"(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. Usefull in simulations with strong velocity contrasts between slow bodies and fast bodies. Only used if oriVerlet=true."))
387+ ((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."))
388+ ((Real,verletDist,((void)"Automatically initialized",-.9),,"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."))
389+ #else
390+ ((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."))
391+ #endif
392 ((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]."))
393 ((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|"))
394 ((int,nBins,5,,"Number of velocity bins for striding. If <=0, bin-less strigin is used (this is however DEPRECATED)."))
395
396=== modified file 'pkg/common/InteractionLoop.cpp'
397--- pkg/common/InteractionLoop.cpp 2011-02-27 14:26:15 +0000
398+++ pkg/common/InteractionLoop.cpp 2011-04-26 10:44:50 +0000
399@@ -103,7 +103,11 @@
400
401 bool swap=false;
402 // IGeomDispatcher
403+ #ifdef ORI_VERLET
404+ if(unlikely(!I->functorCache.geom)){
405+ #else
406 if(unlikely(!I->functorCache.geom || !I->functorCache.phys)){
407+ #endif
408 I->functorCache.geom=geomDispatcher->getFunctor2D(b1_->shape,b2_->shape,swap);
409 // returns NULL ptr if no functor exists; remember that and shortcut
410 if(!I->functorCache.geom) {I->functorCache.geomExists=false; continue; }
411@@ -112,10 +116,10 @@
412 // arguments for the geom functor are in the reverse order (dispatcher would normally call goReverse).
413 // we don't remember the fact that is reverse, so we swap bodies within the interaction
414 // and can call go in all cases
415- if(unlikely(swap)){I->swapOrder(); }
416+ if(unlikely(swap)){I->swapOrder();}
417 // body pointers must be updated, in case we swapped
418- const shared_ptr<Body>& b1=Body::byId(I->getId1(),scene);
419- const shared_ptr<Body>& b2=Body::byId(I->getId2(),scene);
420+ const shared_ptr<Body>& b1=swap?b2_:b1_;
421+ const shared_ptr<Body>& b2=swap?b1_:b2_;
422
423 assert(I->functorCache.geom);
424 bool wasReal=I->isReal();
425
426=== modified file 'pkg/dem/NewtonIntegrator.cpp'
427--- pkg/dem/NewtonIntegrator.cpp 2011-04-22 09:01:59 +0000
428+++ pkg/dem/NewtonIntegrator.cpp 2011-04-26 10:44:50 +0000
429@@ -9,7 +9,6 @@
430 #include<yade/pkg/dem/NewtonIntegrator.hpp>
431 #include<yade/core/Scene.hpp>
432 #include<yade/pkg/dem/Clump.hpp>
433-#include<yade/pkg/common/VelocityBins.hpp>
434 #include<yade/lib/base/Math.hpp>
435
436 YADE_PLUGIN((NewtonIntegrator));
437@@ -66,7 +65,21 @@
438 maxVelocitySq=max(maxVelocitySq,state->vel.squaredNorm());
439 #endif
440 }
441-
442+#ifdef ORI_VERLET
443+void NewtonIntegrator::saveMaximaDisplacement(const shared_ptr<Body>& b){
444+ Vector3r disp=b->state->pos-b->bound->refPos;
445+ Real maxDisp=max(abs(disp[0]),max(abs(disp[1]),abs(disp[2])));
446+ if (!maxDisp || maxDisp<=b->bound->sweepLength) {b->bound->isBounding = (updatingDispFactor>0 && (updatingDispFactor*maxDisp)<b->bound->sweepLength);
447+ maxDisp=0.5;//not 0, else it will be seen as "not updated" by the collider
448+ }
449+ else {b->bound->isBounding = false; maxDisp=2;}
450+ #ifdef YADE_OPENMP
451+ Real& thrMaxVSq=threadMaxVelocitySq[omp_get_thread_num()]; thrMaxVSq=max(thrMaxVSq,maxDisp);
452+ #else
453+ maxVelocitySq=max(maxVelocitySq,maxDisp);
454+ #endif
455+}
456+#endif
457 void NewtonIntegrator::action()
458 {
459 scene->forces.sync();
460@@ -108,21 +121,7 @@
461 if(unlikely(!b || b->isClumpMember())) continue;
462
463 State* state=b->state.get(); const Body::id_t& id=b->getId();
464- Vector3r f=Vector3r::Zero(), m=Vector3r::Zero();
465- // clumps forces
466- if(b->isClump()) {
467- b->shape->cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
468- #ifdef YADE_OPENMP
469- //it is safe here, since onky one thread will read/write
470- scene->forces.getTorqueUnsynced(id)+=m;
471- scene->forces.getForceUnsynced(id)+=f;
472- #else
473- scene->forces.addTorque(id,m);
474- scene->forces.addForce(id, f);
475- #endif
476- }
477- //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)
478- f=scene->forces.getForce(id); m=scene->forces.getTorque(id);
479+ Vector3r f=scene->forces.getForce(id), m=scene->forces.getTorque(id);
480 #ifdef YADE_DEBUG
481 if(isnan(f[0])||isnan(f[1])||isnan(f[2])) throw runtime_error(("NewtonIntegrator: NaN force acting on #"+lexical_cast<string>(id)+".").c_str());
482 if(isnan(m[0])||isnan(m[1])||isnan(m[2])) throw runtime_error(("NewtonIntegrator: NaN torque acting on #"+lexical_cast<string>(id)+".").c_str());
483@@ -134,8 +133,6 @@
484 // in aperiodic boundaries, it is equal to absolute velocity
485 Vector3r fluctVel=isPeriodic?scene->cell->bodyFluctuationVel(b->state->pos,b->state->vel):state->vel;
486
487-
488-
489 // numerical damping & kinetic energy
490 if(unlikely(trackEnergy)) updateEnergy(b,state,fluctVel,f,m);
491
492@@ -144,6 +141,8 @@
493
494 // for particles not totally blocked, compute accelerations; otherwise, the computations would be useless
495 if (state->blockedDOFs!=State::DOF_ALL) {
496+ // forces
497+ if(b->isClump()) b->shape->cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
498 // linear acceleration
499 Vector3r linAccel=computeAccel(f,state->mass,state->blockedDOFs);
500 cundallDamp2nd(dt,f,fluctVel,linAccel);
501@@ -163,7 +162,10 @@
502 leapfrogTranslate(state,id,dt);
503 if(!useAspherical) leapfrogSphericalRotate(state,id,dt);
504 else leapfrogAsphericalRotate(state,id,dt,m);
505-
506+ #ifdef ORI_VERLET
507+ if (oriVerlet) saveMaximaDisplacement(b);
508+ else
509+ #endif
510 saveMaximaVelocity(id,state);
511 // move individual members of the clump, save maxima velocity (for collider stride)
512 if(b->isClump()) Clump::moveMembers(b,scene,this);
513
514=== modified file 'pkg/dem/NewtonIntegrator.hpp'
515--- pkg/dem/NewtonIntegrator.hpp 2011-04-22 09:01:59 +0000
516+++ pkg/dem/NewtonIntegrator.hpp 2011-04-26 10:44:50 +0000
517@@ -11,6 +11,7 @@
518 #include<yade/core/Interaction.hpp>
519 #include<yade/lib/base/Math.hpp>
520 #include<yade/pkg/common/Callbacks.hpp>
521+#include<yade/pkg/common/VelocityBins.hpp>
522 #ifdef YADE_OPENMP
523 #include<omp.h>
524 #endif
525@@ -38,14 +39,19 @@
526
527 // whether the cell has changed from the previous step
528 bool cellChanged;
529+
530 // wether a body has been selected in Qt view
531 bool bodySelected;
532 int homoDeform; // updated from scene at every call; -1 for aperiodic simulations, otherwise equal to scene->cell->homoDeform
533 Matrix3r dVelGrad; // dtto
534
535 public:
536+ Real updatingDispFactor;//(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.
537 // function to save maximum velocity, for the verlet-distance optimization
538 void saveMaximaVelocity(const Body::id_t& id, State* state);
539+ #ifdef ORI_VERLET
540+ void saveMaximaDisplacement(const shared_ptr<Body>& b);
541+ #endif
542 #ifdef YADE_OPENMP
543 vector<Real> threadMaxVelocitySq;
544 #endif
545@@ -56,6 +62,7 @@
546 ((Real,damping,0.2,,"damping coefficient for Cundall's non viscous damping (see [Chareyre2005]_) [-]"))
547 ((Real,maxVelocitySq,NaN,,"store square of max. velocity, for informative purposes; computed again at every step. |yupdate|"))
548 ((bool,exactAsphericalRot,true,,"Enable more exact body rotation integrator for :yref:`aspherical bodies<Body.aspherical>` *only*, using formulation from [Allen1989]_, pg. 89."))
549+ ((bool,oriVerlet,false,,"Compare Verlet distance with displacement instead of velocity (updated by the collider)"))
550 ((Matrix3r,prevVelGrad,Matrix3r::Zero(),,"Store previous velocity gradient (:yref:`Cell::velGrad`) to track acceleration. |yupdate|"))
551 #ifdef YADE_BODY_CALLBACK
552 ((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."))

Subscribers

People subscribed via source and target branches