yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06584
[Branch ~yade-dev/yade/trunk] Rev 2615: - Deeper integration of isDynamic semantic : skip velocity integration, damping, DOFs blocking, a...
------------------------------------------------------------
revno: 2615
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: yade
timestamp: Mon 2010-12-20 16:53:06 +0100
message:
- Deeper integration of isDynamic semantic : skip velocity integration, damping, DOFs blocking, and exact spherical rotations, when !isDynamic. 0-density bodies should be safe.
- Partial DOFs blocking is not implemented yet.
modified:
pkg/dem/NewtonIntegrator.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/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2010-12-19 15:59:14 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2010-12-20 15:53:06 +0000
@@ -11,7 +11,7 @@
#include<yade/pkg/dem/Clump.hpp>
#include<yade/pkg/common/VelocityBins.hpp>
#include<yade/lib/base/Math.hpp>
-
+
YADE_PLUGIN((NewtonIntegrator));
CREATE_LOGGER(NewtonIntegrator);
void NewtonIntegrator::cundallDamp(const Real& dt, const Vector3r& N, const Vector3r& V, Vector3r& A){
@@ -41,7 +41,7 @@
void NewtonIntegrator::handleClumpMemberAngAccel(Scene* scene, const Body::id_t& memberId, State* memberState, State* clumpState){
// angular acceleration from: normal torque + torque generated by the force WRT particle centroid on the clump centroid
const Vector3r& m=scene->forces.getTorque(memberId); const Vector3r& f=scene->forces.getForce(memberId);
- Vector3r diffClumpAngularAccel=m.cwise()/clumpState->inertia+(memberState->pos-clumpState->pos).cross(f).cwise()/clumpState->inertia;
+ Vector3r diffClumpAngularAccel=m.cwise()/clumpState->inertia+(memberState->pos-clumpState->pos).cross(f).cwise()/clumpState->inertia;
// damp increment of accels on the clump, using velocities of the clump MEMBER
cundallDamp(scene->dt,m,memberState->angVel,diffClumpAngularAccel);
clumpState->angAccel+=diffClumpAngularAccel;
@@ -99,7 +99,7 @@
State* state=b->state.get();
const Body::id_t& id=b->getId();
// clump members are non-dynamic; we only get their velocities here
- if(unlikely(b->isClumpMember())){
+ if(unlikely(/*!b->isDynamic() ||*/ b->isClumpMember())){
saveMaximaVelocity(scene,id,state);
continue;
}
@@ -137,14 +137,16 @@
if (likely(b->isStandalone())){
// translate equation
- state->accel=f/state->mass;
- cundallDamp(dt,f,fluctVel,state->accel);
+ if(b->isDynamic()){
+ state->accel=f/state->mass;
+ cundallDamp(dt,f,fluctVel,state->accel);}
leapfrogTranslate(scene,state,id,dt);
// rotate equation
// exactAsphericalRot is disabled or the body is spherical
- if (likely(!exactAsphericalRot || !b->isAspherical())){
- state->angAccel=m.cwise()/state->inertia;
- cundallDamp(dt,m,state->angVel,state->angAccel);
+ if (likely(!exactAsphericalRot || !b->isAspherical() || !b->isDynamic())){
+ if(b->isDynamic()) {
+ state->angAccel=m.cwise()/state->inertia;
+ cundallDamp(dt,m,state->angVel,state->angAccel);}
leapfrogSphericalRotate(scene,state,id,dt);
} else { // exactAsphericalRot enabled & aspherical body
// no damping in this case
@@ -173,9 +175,10 @@
leapfrogTranslate(scene,state,id,dt);
leapfrogAsphericalRotate(scene,state,id,dt,M);
} else { // exactAsphericalRot disabled or clump is spherical
- Vector3r dAngAccel=M.cwise()/state->inertia;
- cundallDamp(dt,M,state->angVel,dAngAccel);
- state->angAccel+=dAngAccel;
+ if(b->isDynamic()){
+ Vector3r dAngAccel=M.cwise()/state->inertia;
+ cundallDamp(dt,M,state->angVel,dAngAccel);
+ state->angAccel+=dAngAccel;}
FOREACH(Clump::MemberMap::value_type mm, static_cast<Clump*>(b->shape.get())->members){
const Body::id_t& memberId=mm.first;
const shared_ptr<Body>& member=Body::byId(memberId,scene); assert(member->isClumpMember());
@@ -207,38 +210,39 @@
inline void NewtonIntegrator::leapfrogTranslate(Scene* scene, State* state, const Body::id_t& id, const Real& dt)
{
- blockTranslateDOFs(state->blockedDOFs, state->accel);
-
+ if(state->blockedDOFs!=State::DOF_ALL) blockTranslateDOFs(state->blockedDOFs, state->accel);
+
if (scene->forces.getMoveRotUsed()) state->pos+=scene->forces.getMove(id);
if (homoDeform==Cell::HOMO_VEL || homoDeform==Cell::HOMO_VEL_2ND) {
// update velocity reflecting changes in the macroscopic velocity field, making the problem homothetic.
//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).
- //This is the convective term, appearing in the time derivation of Cundall/Thornton expression (dx/dt=velGrad*pos -> d²x/dt²=dvelGrad/dt+velGrad*vel), negligible in many cases but not for high speed large deformations (gaz or turbulent flow). Emulating Cundall is an option, I don't especially recommend it. I know homothetic 1 and 2 expressions tend to identical values in the limit of dense quasi-static situations. They can give slitghly different results in other cases, and I'm still not sure which one should be considered better, if any (Cundall formula is in contradiction with molecular dynamics litterature).
+ //This is the convective term, appearing in the time derivation of Cundall/Thornton expression (dx/dt=velGrad*pos -> d²x/dt²=dvelGrad/dt+velGrad*vel), negligible in many cases but not for high speed large deformations (gaz or turbulent flow).
if (homoDeform==Cell::HOMO_VEL_2ND) state->vel+=scene->cell->prevVelGrad*state->vel*dt;
-
+
//In all cases, reflect macroscopic (periodic cell) acceleration in the velocity. This is the dominant term in the update in most cases
Vector3r dVel=dVelGrad*state->pos;
state->vel+=dVel;
} else if (homoDeform==Cell::HOMO_POS){
state->pos+=scene->cell->velGrad*state->pos*dt;
}
- state->vel+=dt*state->accel;
+ if(state->blockedDOFs!=State::DOF_ALL) state->vel+=dt*state->accel;
state->pos += state->vel*dt;
}
inline void NewtonIntegrator::leapfrogSphericalRotate(Scene* scene, State* state, const Body::id_t& id, const Real& dt )
{
- blockRotateDOFs(state->blockedDOFs, state->angAccel);
- state->angVel+=dt*state->angAccel;
+ if(state->blockedDOFs!=State::DOF_ALL){
+ blockRotateDOFs(state->blockedDOFs, state->angAccel);
+ state->angVel+=dt*state->angAccel;}
Vector3r axis = state->angVel;
-
+
if (axis!=Vector3r::Zero()) { //If we have an angular velocity, we make a rotation
Real angle=axis.norm(); axis/=angle;
Quaternionr q(AngleAxisr(angle*dt,axis));
state->ori = q*state->ori;
}
-
+
if(scene->forces.getMoveRotUsed() && scene->forces.getRot(id)!=Vector3r::Zero()) {
Vector3r r(scene->forces.getRot(id));
Real norm=r.norm(); r/=norm;
@@ -269,7 +273,7 @@
Quaternionr q(AngleAxisr(norm,r));
state->ori=q*state->ori;
}
- state->ori.normalize();
+ state->ori.normalize();
}
// http://www.euclideanspace.com/physics/kinematics/angularvelocity/QuaternionDifferentiation2.pdf
Quaternionr NewtonIntegrator::DotQ(const Vector3r& angVel, const Quaternionr& Q){