>> Now, more annoying : type this once the beam is stabilized :
>> O.bodies[0].state.angVel=Vector3(0,0,0.01)
Actually, this trick doesn't work any longer after r2607 ("Remove
Body::flags/FLAG_DYNAMIC, setDynamic now merely sets
state->blockedDOFs=State::DOF_ALL etc.") : Newton will totally skip non
dynamic bodies. I've been really scared for a few minutes looking at
absurd beam motion, but finally the fix in r2608 is correct.
I have fixs in Newton to enables manual velocity setting again (together
with optimized non dynamic bodies handling), but it is only partialy
implemented (not sure about clumps and nonspherical integration), and
more annoying, it fails at regression (below). The diff is attached, I
won't commit.
Bruno
======================================================================
FAIL: testTranslationRotationEngines (yade.TestEngines)
----------------------------------------------------------------------
Traceback (most recent call last):
File
"/home/3S-LAB/bchareyre/yade/yade-test-kdev/lib/yade-true-mono/py/yade/tests/omega.py",
line 85, in testTranslationRotationEngines
self.assertTrue(int(O.bodies[id_nodyn_transl].state.pos[0]) ==
O.iter)
#Check translation of nondynamic bodies
AssertionError
>> O.bodies[0].state.angVel=Vector3(0,0,0.01)
>> In both cases you will see growing errors in the form of a weird deformed
>> shape after
>> enough 360° turns (earlier with ScGeom, I confess).
>
> Fixed in 2608. McNamara win a battle.
>
> Bruno
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : [email protected]
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
>
--
_______________
Bruno Chareyre
Associate Professor
ENSE³ - Grenoble INP
Lab. 3SR
BP 53 - 38041, Grenoble cedex 9 - France
Tél : +33 4 56 52 86 21
Fax : +33 4 76 82 70 43
________________
=== modified file 'pkg/dem/NewtonIntegrator.cpp'
--- pkg/dem/NewtonIntegrator.cpp 2010-12-05 17:10:06 +0000
+++ pkg/dem/NewtonIntegrator.cpp 2010-12-10 20:09:44 +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->isDynamic() || 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(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,8 +210,8 @@
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.
@@ -216,29 +219,30 @@
//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).
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){
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : [email protected]
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
