← Back to team overview

yade-users team mailing list archive

  1. yade-users team
  2. Mailing list archive
  3. Message #01045

Re: How to reset a body's velocity?

  Thread Previous • Date Previous • • Thread Next 
sega wrote:
Actually, the fact that prevVelocities are inside the integrator could be considered a flaw in design... 

Hmm... As me seems, prevVelocities is specific for leapfrog method.
For another integrator may be another implementation (Gear predictor-corrector as example...)


Ok...

My opinion today is :

1- Leapfrog integrator does NOT need a separate variable for
previousVelocity, since at each timestep :

 vel = vel + acceleration*dt; //which is equivalent with NewVel =
PreviousVel + acceleration*dt

2- The current LeapFrog equations in Yade are not consistent, because
for translation we have :
/rb->se3.position += prevVelocities[id]*dt;/
and for rotation we have :
/rb->se3.orientation = q*rb->se3.orientation/; //with q defined _based
on rb->angularVelocity_

This results in position X(t+dt) being computed as
X(t)+dX/dt(t+dt/2)*dt, while orientation O(t+dt)=O(t)+dO/dt(t)*dt.
This means that translation is integrated with a 2nd order scheme (correct) while 
orientation uses a 1st order formulation!!!
Perhaps not really wrong, but surely not optimal.


3- In the contact law, relative displacement is defined based on
body->velocity. Consequence : it does not correspond to the actual
displacement during dt (remember that the displacement is
prevVelocity*dt). This is wrong (or... lets say... one half of order
scheme...).

4- I removed previousVel/previousAngularVel, to get back consistency for
translation/rotation and for body motion/contact kinematics.
Result :
- the scheme is a lot more stable, for the first time in Yade I get
unbalanced force around 10e-7
- and this stability remains true even with dt 6 times larger than before.

Thank you Sega for pointing that out.
I will think about that a bit more before commiting anything, also
please let me know your opinion about this.
Sega : one side effect will be that body->velocity=0 will be enough to
reset a body velocity. :)

Bruno

p.s. I'm sending a modified NewtonsDampedLaw attached to this mail,
previousVel is still in the code but its overwritten with
body->velocity. Could you guys please test it and confirm that it works
well in your simulations?









/*************************************************************************
 Copyright (C) 2008 by Bruno Chareyre		                         *
*  bruno.chareyre@xxxxxxxxxxx      					 *
*                                                                        *
*  This program is free software; it is licensed under the terms of the  *
*  GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/

#include "NewtonsDampedLaw.hpp"
#include<yade/core/MetaBody.hpp>
#include <yade/pkg-common/RigidBodyParameters.hpp>
#include <yade/pkg-common/Momentum.hpp>
#include <yade/pkg-common/Force.hpp>
#include<yade/lib-base/yadeWm3Extra.hpp>


void NewtonsDampedLaw::registerAttributes()
{
	//DeusExMachina::registerAttributes();
	REGISTER_ATTRIBUTE(damping);
}

NewtonsDampedLaw::NewtonsDampedLaw()
{
	damping = 0.2;
	forceClassIndex = (new Force)->getClassIndex();
	momentumClassIndex = (new Momentum)->getClassIndex();
	prevSize = 0;
	
}

void NewtonsDampedLaw::applyCondition ( Body * body )
{


	MetaBody * ncb = static_cast<MetaBody*> ( body );
	shared_ptr<BodyContainer>& bodies = ncb->bodies;
	Real dt = Omega::instance().getTimeStep();
	BodyContainer::iterator bi    = bodies->begin();
	BodyContainer::iterator biEnd = bodies->end();
	for ( ; bi!=biEnd ; ++bi )
	{

		const shared_ptr<Body>& b = *bi;

		if (!b->isDynamic) continue;
		
		RigidBodyParameters * rb = YADE_CAST<RigidBodyParameters*> ( b->physicalParameters.get() );
		unsigned int id = b->getId();
		Vector3r& m = ( static_cast<Momentum*> ( ncb->physicalActions->find ( id, momentumClassIndex ).get() ) )->momentum;
		Vector3r& f = ( static_cast<Force*> ( ncb->physicalActions->find ( id, forceClassIndex ).get() ) )->force;

		//Orientation and Position Integrators :
		if ( prevSize <=id )
		{
			prevSize = id+1;
			prevAngularVelocities.resize ( prevSize );
			prevVelocities.resize ( prevSize );
			firsts.resize ( prevSize,true );
		}

		if (!b->isClump()) {
			//Damping moments
			for ( int i=0; i<3; ++i )
			{
				m[i] *= 1 - damping*(firsts[id] ? 0 : Mathr::Sign ( m[i]*(prevAngularVelocities[id][i]+m[i]/rb->inertia[i]*0.5*dt)));
			}
			//Damping force :
			for ( int i=0; i<3; ++i )
			{
				f[i] *= 1 - damping*(firsts[id] ? 0 : Mathr::Sign ( f[i]*(prevVelocities[id][i] + ( ( Real ) 0.5 ) *dt*f[i]/rb->mass)));
			}
		}
		//Newtons mometum law :
		if ( b->isStandalone() ) rb->angularAcceleration=diagDiv ( m,rb->inertia );
		else if ( b->isClump() ) rb->angularAcceleration+=diagDiv ( m,rb->inertia );
		else
		{ // isClumpMember()
			const shared_ptr<Body>& clump ( Body::byId ( b->clumpId ) );
			RigidBodyParameters* clumpRBP=YADE_CAST<RigidBodyParameters*> ( clump->physicalParameters.get() );
			/* angularAcceleration is reset by ClumpMemberMover engine */
			clumpRBP->angularAcceleration+=diagDiv ( m,clumpRBP->inertia );
		}

		// Newtons force law
		if ( b->isStandalone() ) rb->acceleration=f/rb->mass;
		else if ( b->isClump() )
		{
			// accel for clump reset in Clump::moveMembers, called by ClumpMemberMover engine
			rb->acceleration+=f/rb->mass;
		}
		else
		{ // force applied to a clump member is applied to clump itself
			const shared_ptr<Body>& clump ( Body::byId ( b->clumpId ) );
			RigidBodyParameters* clumpRBP=YADE_CAST<RigidBodyParameters*> ( clump->physicalParameters.get() );
			// accels reset by Clump::moveMembers in last iteration
			clumpRBP->acceleration+=f/clumpRBP->mass;
			clumpRBP->angularAcceleration+=diagDiv ( ( rb->se3.position-clumpRBP->se3.position ).Cross ( f ),clumpRBP->inertia ); //acceleration from torque generated by the force WRT particle centroid on the clump centroid
		}

		

		

		if ( !firsts[id] )
		{
			rb->angularVelocity = prevAngularVelocities[id]+rb->angularAcceleration*0.5*dt;
			rb->velocity = prevVelocities[id]+ ( ( Real ) 0.5 ) *dt*rb->acceleration;
		}

		prevAngularVelocities[id] = rb->angularVelocity+ ( ( Real ) ( 0.5*dt ) ) *rb->angularAcceleration;
		prevVelocities[id] = rb->velocity+ ( ( Real ) 0.5 ) *dt*rb->acceleration;
		rb->velocity = prevVelocities[id];
		rb->angularVelocity = prevAngularVelocities[id];

		//Vector3r axis2 = rb->angularVelocity;
		Vector3r axis = prevAngularVelocities[id];
		//cerr << "axis = "<< axis << " axis2 =" << axis2 << endl;
		Real angle = axis.Normalize();
		Quaternionr q;
		q.FromAxisAngle ( axis,angle*dt );
		rb->se3.orientation = q*rb->se3.orientation;
		//cerr << "rb->se3.orientation= " << rb->se3.orientation[0] << " " << rb->se3.orientation[1] << " "<< endl;
		rb->se3.orientation.Normalize();


		rb->se3.position += prevVelocities[id]*dt;

		firsts[id] = false;


	}
}

YADE_PLUGIN();


_______________________________________________
Yade-users mailing list
Yade-users@xxxxxxxxxxxxxxxx
https://lists.berlios.de/mailman/listinfo/yade-users

Follow ups

References

  Thread Previous • Date Previous • • Thread Next