yade-dev team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

>> 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     : yade-dev@xxxxxxxxxxxxxxxxxxx
> 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){