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