yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #08259
[Branch ~yade-dev/yade/trunk] Rev 3007: - Collider: fix the probe function (account for current body's position + typo in one comparison)
------------------------------------------------------------
revno: 3007
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Mon 2012-01-30 11:42:13 +0100
message:
- Collider: fix the probe function (account for current body's position + typo in one comparison)
- Newton: remove the useless "force" from damping2nd parameters, don't apply meanfield vel when not periodic (only diff is it avoids a multiplication of null matrix and vector), make gravity attribute work correctly
- GravityEngine: add LOG_ERROR on deprecation (LOG_WARNING are not displayed in prebuilt package right?), better use Newton::gravity to save cpu time.
- Clump.hpp: fix the template function that was using the old saveMaximaVelocity, use the new saveMaximaDisplacement instead
- PeriIsoCompressor: never pause() in some c++ code, since it resuults in stuck GUI without a chance to run more timesteps
modified:
pkg/common/GravityEngines.cpp
pkg/common/GravityEngines.hpp
pkg/common/InsertionSortCollider.cpp
pkg/common/InsertionSortCollider.hpp
pkg/dem/Clump.hpp
pkg/dem/NewtonIntegrator.cpp
pkg/dem/NewtonIntegrator.hpp
pkg/dem/PeriIsoCompressor.cpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/common/GravityEngines.cpp'
--- pkg/common/GravityEngines.cpp 2011-11-30 17:39:33 +0000
+++ pkg/common/GravityEngines.cpp 2012-01-30 10:42:13 +0000
@@ -13,8 +13,10 @@
#include<boost/regex.hpp>
YADE_PLUGIN((GravityEngine)(CentralGravityEngine)(AxialGravityEngine)(HdapsGravityEngine));
+CREATE_LOGGER(GravityEngine);
void GravityEngine::action(){
+ if (warnOnce) {warnOnce=false; LOG_ERROR("GravityEngine is deprecated, consider using Newton::gravity instead (unless gravitational energy has to be tracked - not implemented with the newton attribute).")}
const bool trackEnergy(unlikely(scene->trackEnergy));
const Real dt(scene->dt);
YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>& b, scene->bodies){
=== modified file 'pkg/common/GravityEngines.hpp'
--- pkg/common/GravityEngines.hpp 2011-03-01 10:31:54 +0000
+++ pkg/common/GravityEngines.hpp 2012-01-30 10:42:13 +0000
@@ -13,8 +13,10 @@
((Vector3r,gravity,Vector3r::Zero(),,"Acceleration [kgms⁻²]"))
((int,gravPotIx,-1,(Attr::noSave|Attr::hidden),"Index for gravPot energy"))
((int,mask,0,,"If mask defined, only bodies with corresponding groupMask will be affected by this engine. If 0, all bodies will be affected."))
+ ((bool,warnOnce,true,,"For deprecation warning once."))
,/*ctor*/,/*py*/
);
+ DECLARE_LOGGER;
};
REGISTER_SERIALIZABLE(GravityEngine);
=== modified file 'pkg/common/InsertionSortCollider.cpp'
--- pkg/common/InsertionSortCollider.cpp 2012-01-26 11:19:05 +0000
+++ pkg/common/InsertionSortCollider.cpp 2012-01-30 10:42:13 +0000
@@ -25,15 +25,6 @@
///fast
if (spatialOverlap(id1,id2) && Collider::mayCollide(Body::byId(id1,scene).get(),Body::byId(id2,scene).get()) && !interactions->found(id1,id2))
interactions->insert(shared_ptr<Interaction>(new Interaction(id1,id2)));
-
- ///slow
-// //existing interaction?
-// if (interactions->found(id1,id2)) {
-// //if existing interaction is virtual and bounds don't overlap, remove it
-// if (!interactions->find(id1,id2)->isReal() && !spatialOverlap(id1,id2)){ interactions->erase(id1,id2); return;}}
-// //if it doesn't exist and bounds overlap, create a virtual interaction
-// else if (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){
@@ -74,13 +65,14 @@
int offset = 3*it->id;
const shared_ptr<Body>& b=Body::byId(it->id,scene);
if(unlikely(!b)) continue;
- const Real& sweepLength = b->bound->sweepLength;
- if (!(maxima[offset]-sweepLength < bv.min[0] ||
- minima[offset]+sweepLength > bv.min[0] ||
- minima[offset+1]+sweepLength > bv.max[1] ||
- maxima[offset+1]-sweepLength < bv.min[1] ||
- minima[offset+2]+sweepLength > bv.max[2] ||
- maxima[offset+2]-sweepLength < bv.min[2] ))
+ const Real& sweepLength = b->bound->sweepLength;
+ Vector3r disp = b->state->pos - b->bound->refPos;
+ if (!(maxima[offset]-sweepLength+disp[0] < bv.min[0] ||
+ minima[offset]+sweepLength+disp[0] > bv.max[0] ||
+ minima[offset+1]+sweepLength+disp[1] > bv.max[1] ||
+ maxima[offset+1]-sweepLength+disp[1] < bv.min[1] ||
+ minima[offset+2]+sweepLength+disp[2] > bv.max[2] ||
+ maxima[offset+2]-sweepLength+disp[2] < bv.min[2] ))
{
ret.push_back(it->id);
}
@@ -98,9 +90,6 @@
if(fastestBodyMaxDist>=1 || fastestBodyMaxDist==0) return true;
if((size_t)BB[0].size!=2*scene->bodies->size()) return true;
if(scene->interactions->dirty) return true;
- // we wouldn't run in this step; in that case, just delete pending interactions
- // this is done in ::action normally, but it would make the call counters not reflect the stride
-// scene->interactions->erasePending(*this,scene);
return false;
}
@@ -156,6 +145,7 @@
if(!s) continue;
minR=min(s->radius,minR);
}
+ if (isinf(minR)) LOG_ERROR("verletDist is set to 0 because no spheres were found. It will result in suboptimal performances, consider setting a positive verletDist in your script.");
// if no spheres, disable stride
verletDist=isinf(minR) ? 0 : abs(verletDist)*minR;
}
=== modified file 'pkg/common/InsertionSortCollider.hpp'
--- pkg/common/InsertionSortCollider.hpp 2012-01-23 14:43:54 +0000
+++ pkg/common/InsertionSortCollider.hpp 2012-01-30 10:42:13 +0000
@@ -209,10 +209,10 @@
((int,sortAxis,0,,"Axis for the initial contact detection."))
((bool,allowBiggerThanPeriod,false,,"If true, tests on bodies sizes will be disabled, and the simulation will run normaly even if bodies larger than period are found. It can be usefull when the periodic problem include e.g. a floor modelized with wall/box/facet.\nBe sure you know what you are doing if you touch this flag. The result is undefined if one large body moves out of the (0,0,0) period."))
((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.)"))
- ((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."))
+ ((int,targetInterv,50,,"(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."))
((Real,updatingDispFactor,-1,,"(experimental) Displacement factor used to trigger bound update: the bound is updated only if updatingDispFactor*disp>sweepDist when >0, else all bounds are updated."))
((Real,verletDist,((void)"Automatically initialized",-.5),,"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. The actual length added to one bound can be only a fraction of verletDist when :yref:`InsertionSortCollider::targetInterv` is >0."))
- ((Real,minSweepDistFactor,0.2,,"Minimal distance by which enlarge all bounding boxes; superseeds computed value of verletDist when lower that (minSweepDistFactor x verletDist)."))
+ ((Real,minSweepDistFactor,0.1,,"Minimal distance by which enlarge all bounding boxes; superseeds computed value of verletDist when lower that (minSweepDistFactor x verletDist)."))
((Real,fastestBodyMaxDist,-1,,"Normalized maximum displacement of the fastest body since last run; if >= 1, we could get out of bboxes and will trigger full run. |yupdate|"))
((int,numReinit,0,Attr::readonly,"Cummulative number of bound array re-initialization."))
((Real,useless,,,"for compatibility of scripts defining the old collider's attributes - see deprecated attributes"))
=== modified file 'pkg/dem/Clump.hpp'
--- pkg/dem/Clump.hpp 2011-01-10 08:35:00 +0000
+++ pkg/dem/Clump.hpp 2012-01-30 10:42:13 +0000
@@ -70,14 +70,15 @@
const shared_ptr<State>& clumpState=clumpBody->state;
FOREACH(MemberMap::value_type& B, clump->members){
// B.first is Body::id_t, B.second is local Se3r of that body in the clump
- const shared_ptr<State>& subState=Body::byId(B.first,scene)->state; const Vector3r& subPos(B.second.position); const Quaternionr& subOri(B.second.orientation);
+ const shared_ptr<Body>& b = Body::byId(B.first,scene);
+ const shared_ptr<State>& subState=b->state; const Vector3r& subPos(B.second.position); const Quaternionr& subOri(B.second.orientation);
// position update
subState->pos=clumpState->pos+clumpState->ori*subPos;
subState->ori=clumpState->ori*subOri;
// velocity update
subState->vel=clumpState->vel+clumpState->angVel.cross(subState->pos-clumpState->pos);
subState->angVel=clumpState->angVel;
- if(likely(integrator)) integrator->saveMaximaVelocity(B.first,subState.get());
+ if(likely(integrator)) integrator->saveMaximaDisplacement(b);
}
}
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2012-01-24 03:29:29 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2012-01-30 10:42:13 +0000
@@ -19,14 +19,14 @@
for(int i=0; i<3; i++) force[i]*=1-damping*Mathr::Sign(force[i]*vel[i]);
}
// 2nd order numerical damping
-void NewtonIntegrator::cundallDamp2nd(const Real& dt, const Vector3r& force, const Vector3r& vel, Vector3r& accel){
- for(int i=0; i<3; i++) accel[i]*= 1 - damping*Mathr::Sign ( force[i]*(vel[i] + 0.5*dt*accel[i]) );
+void NewtonIntegrator::cundallDamp2nd(const Real& dt, const Vector3r& vel, Vector3r& accel){
+ for(int i=0; i<3; i++) accel[i]*= 1 - damping*Mathr::Sign ( accel[i]*(vel[i] + 0.5*dt*accel[i]) );
}
Vector3r NewtonIntegrator::computeAccel(const Vector3r& force, const Real& mass, int blockedDOFs){
- if(likely(blockedDOFs==0)) return force/mass;
+ if(likely(blockedDOFs==0)) return (force/mass + gravity);
Vector3r ret(Vector3r::Zero());
- for(int i=0; i<3; i++) if(!(blockedDOFs & State::axisDOF(i,false))) ret[i]+=force[i]/mass;
+ for(int i=0; i<3; i++) if(!(blockedDOFs & State::axisDOF(i,false))) ret[i]+=force[i]/mass+gravity[i];
return ret;
}
Vector3r NewtonIntegrator::computeAngAccel(const Vector3r& torque, const Vector3r& inertia, int blockedDOFs){
@@ -119,7 +119,7 @@
if(unlikely(b->isClumpMember())) continue;
State* state=b->state.get(); const Body::id_t& id=b->getId();
- Vector3r f=gravity*state->mass, m=Vector3r::Zero();
+ Vector3r f=Vector3r::Zero(), m=Vector3r::Zero();
// clumps forces
if(b->isClump()) {
b->shape->cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
@@ -145,19 +145,17 @@
// in aperiodic boundaries, it is equal to absolute velocity
Vector3r fluctVel=isPeriodic?scene->cell->bodyFluctuationVel(b->state->pos,b->state->vel,prevVelGrad):state->vel;
-
-
// numerical damping & kinetic energy
if(unlikely(trackEnergy)) updateEnergy(b,state,fluctVel,f,m);
// whether to use aspherical rotation integration for this body; for no accelerations, spherical integrator is "exact" (and faster)
- bool useAspherical=(b->isAspherical() && exactAsphericalRot && state->blockedDOFs!=State::DOF_ALL);
+ bool useAspherical=(exactAsphericalRot && b->isAspherical() && state->blockedDOFs!=State::DOF_ALL);
// for particles not totally blocked, compute accelerations; otherwise, the computations would be useless
if (state->blockedDOFs!=State::DOF_ALL) {
// linear acceleration
Vector3r linAccel=computeAccel(f,state->mass,state->blockedDOFs);
- if(state->isDamped) cundallDamp2nd(dt,f,fluctVel,linAccel);
+ if(state->isDamped) cundallDamp2nd(dt,fluctVel,linAccel);
//This is the convective term, appearing in the time derivation of Cundall/Thornton expression (dx/dt=velGrad*pos -> d²x/dt²=dvelGrad/dt*pos+velGrad*vel), negligible in many cases but not for high speed large deformations (gaz or turbulent flow).
linAccel+=prevVelGrad*state->vel;
//finally update velocity
@@ -165,7 +163,7 @@
// angular acceleration
if(!useAspherical){ // uses angular velocity
Vector3r angAccel=computeAngAccel(m,state->inertia,state->blockedDOFs);
- if(state->isDamped) cundallDamp2nd(dt,m,state->angVel,angAccel);
+ if(state->isDamped) cundallDamp2nd(dt,state->angVel,angAccel);
state->angVel+=dt*angAccel;
} else { // uses torque
for(int i=0; i<3; i++) if(state->blockedDOFs & State::axisDOF(i,true)) m[i]=0; // block DOFs here
@@ -201,8 +199,7 @@
//NOTE : if the velocity is updated before moving the body, it means the current velGrad (i.e. before integration in cell->integrateAndUpdate) will be effective for the current time-step. Is it correct? If not, this velocity update can be moved just after "state->pos += state->vel*dt", meaning the current velocity impulse will be applied at next iteration, after the contact law. (All this assuming the ordering is resetForces->integrateAndUpdate->contactLaw->PeriCompressor->NewtonsLaw. Any other might fool us.)
//NOTE : dVel defined without wraping the coordinates means bodies out of the (0,0,0) period can move realy fast. It has to be compensated properly in the definition of relative velocities (see Ig2 functors and contact laws).
//Reflect mean-field (periodic cell) acceleration in the velocity
- Vector3r dVel=dVelGrad*state->pos;
- state->vel+=dVel;
+ if(scene->isPeriodic) {Vector3r dVel=dVelGrad*state->pos; state->vel+=dVel;}
if (!bodySelected || scene->selectedBody!=id) state->pos+=state->vel*dt;
}
=== modified file 'pkg/dem/NewtonIntegrator.hpp'
--- pkg/dem/NewtonIntegrator.hpp 2012-01-23 14:43:54 +0000
+++ pkg/dem/NewtonIntegrator.hpp 2012-01-30 10:42:13 +0000
@@ -22,7 +22,7 @@
class NewtonIntegrator : public GlobalEngine{
inline void cundallDamp1st(Vector3r& force, const Vector3r& vel);
- inline void cundallDamp2nd(const Real& dt, const Vector3r& force, const Vector3r& vel, Vector3r& accel);
+ inline void cundallDamp2nd(const Real& dt, const Vector3r& vel, Vector3r& accel);
inline void leapfrogTranslate(State*, const Body::id_t& id, const Real& dt); // leap-frog translate
inline void leapfrogSphericalRotate(State*, const Body::id_t& id, const Real& dt); // leap-frog rotate of spherical body
inline void leapfrogAsphericalRotate(State*, const Body::id_t& id, const Real& dt, const Vector3r& M); // leap-frog rotate of aspherical body
=== modified file 'pkg/dem/PeriIsoCompressor.cpp'
--- pkg/dem/PeriIsoCompressor.cpp 2011-02-09 09:46:15 +0000
+++ pkg/dem/PeriIsoCompressor.cpp 2012-01-30 10:42:13 +0000
@@ -228,7 +228,7 @@
if (!doneHook.empty()){
LOG_DEBUG ( "Running doneHook: "<<doneHook );
pyRunString(doneHook);}
- else { Omega::instance().pause(); }
+// else { Omega::instance().pause(); }
}
}
}
