yade-dev team mailing list archive

[Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>]

Ok, I understand.
I like your diff more than the previous one. :)
Still one suggestion: this additional parameter in the damping functions
is useless:

+					cundallDamp2nd(dt,m,state->angVel,angAccel,dampcoeff);

In your code, dampCoeff will be either Newton::damping or zero.
So, it would be better to just write this, avoiding the function call
and additional parameter:

if (state->isDamped) cundallDamp2nd(dt,m,state->angVel,angAccel)

Then I agree to commit.

Bruno




On 23/01/12 11:48, Klaus Thoeni wrote:
> Hi Chiara,
>
> thanks for your comment. I am not after viscous damping either. But I think 
> the current implementation does exactly what I want, it's the same as the PFC 
> option b_damp(). So I can deactivate the damping coefficient for individual 
> bodies.
>
> I can commit the code if anyone is interested , after the diff has been 
> approved ;-)
>
> Klaus
>
> On Mon, 23 Jan 2012 08:44:02 PM Chiara wrote:
>> On 23/01/12 00:12, Klaus Thoeni wrote:
>>> Hi Bruno, hi Anton,
>>>
>>> I think density scaling is not what I am looking for since I am
>>> interested in real dynamic simulations. I need different damping
>>> parameters in NewtonIntegrator for the block and the mesh. I have to
>>> consider free flight under gravity so damping=0 (and I am interested in
>>> the real flight time). For the mesh I have to consider additional energy
>>> absorption which is not considered in my model via damping.
>> If you are not after any density scaling, why damping in
>> NewtonIntegrator should be acceptable in a realistic dynamic situation?
>> Would not be easier/more meaningful for you to implement damping inside
>> your contact model? Sorry, just curious.
>> Chiara
>>
>>> It is similar to b_damp() in PFC which is used to remove damping for
>>> certain particle. I implemented it the same way now. The state has a
>>> member which is true by default which means that damping is used. It can
>>> be set to false and damping=0 will be used for this particle. So
>>> basically damping in NewtonIntegrator can be activated and deactivated
>>> for individual particles.
>>>
>>> @Anton: It has to be checked for all particles.
>>>
>>> Here is the diff of my latest implementation, tell me what you think:
>>>
>>> === modified file core/State.hpp
>>> --- core/State.hpp	2011-02-14 08:05:09 +0000
>>> +++ core/State.hpp	2012-01-22 23:56:31 +0000
>>> @@ -59,7 +59,8 @@
>>>
>>>   		((Vector3r,inertia,Vector3r::Zero(),,"Inertia of associated body, 
> in
>>> local coordinate system."))
>>>
>>>   		((Vector3r,refPos,Vector3r::Zero(),,"Reference position"))
>>>   		((Quaternionr,refOri,Quaternionr::Identity(),,"Reference
>>>   		orientation"))
>>>
>>> -		((unsigned,blockedDOFs,,,"[Will be overridden]")),
>>> +		((unsigned,blockedDOFs,,,"[Will be overridden]"))
>>> +		((bool,isDamped,true,,"Damping in :yref:`Newtonintegrator` can be
>>> deactivated for individual particles by setting this variable to FALSE.
>>> E.g. damping is inappropriate for particles in free flight under gravity
>>> but it might still be applicable to other particles in the same
>>> simulation.")),
>>>
>>>   		/* additional initializers */
>>>   		
>>>   			((pos,se3.position))
>>>   			((ori,se3.orientation)),
>>>
>>> === modified file pkg/dem/NewtonIntegrator.cpp
>>> --- pkg/dem/NewtonIntegrator.cpp	2011-11-30 17:39:33 +0000
>>> +++ pkg/dem/NewtonIntegrator.cpp	2012-01-22 23:56:31 +0000
>>> @@ -11,17 +11,18 @@
>>>
>>>   #include<yade/pkg/dem/Clump.hpp>
>>>   #include<yade/pkg/common/VelocityBins.hpp>
>>>   #include<yade/lib/base/Math.hpp>
>>>   
>>>   YADE_PLUGIN((NewtonIntegrator));
>>>   CREATE_LOGGER(NewtonIntegrator);
>>>   
>>>   // 1st order numerical damping
>>>
>>> -void NewtonIntegrator::cundallDamp1st(Vector3r&  force, const Vector3r& 
>>> vel){ -	for(int i=0; i<3; i++)
>>> force[i]*=1-damping*Mathr::Sign(force[i]*vel[i]); +void
>>> NewtonIntegrator::cundallDamp1st(Vector3r&  force, const Vector3r&  vel,
>>> const Real&  dampcoeff){
>>> +	for(int i=0; i<3; i++)
>>> force[i]*=1-dampcoeff*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& 
>>> force, const Vector3r&  vel, Vector3r&  accel, const Real&  dampcoeff){
>>> +	for(int i=0; i<3; i++) accel[i]*= 1 - dampcoeff*Mathr::Sign (
>>> force[i]*(vel[i] + 0.5*dt*accel[i]) );
>>>
>>>   }
>>>   
>>>   Vector3r NewtonIntegrator::computeAccel(const Vector3r&  force, const
>>>   Real&
>>>
>>> mass, int blockedDOFs){
>>> @@ -39,11 +40,13 @@
>>>
>>>   void NewtonIntegrator::updateEnergy(const shared_ptr<Body>&  b, const
>>>   State*
>>>
>>> state, const Vector3r&  fluctVel, const Vector3r&  f, const Vector3r& 
>>> m){
>>>
>>>   	assert(b->isStandalone() || b->isClump());
>>>
>>> -	// always positive dissipation, by-component: |F_i|*|v_i|*damping*dt 
> (|
>>> T_i|*|ω_i|*damping*dt for rotations)
>>> -	if(damping!=0.){
>>> -		scene->energy-
>>>
>>>> add(fluctVel.cwise().abs().dot(f.cwise().abs())*damping*scene-
>>>> dt,"nonviscDamp",nonviscDampIx,/*non-incremental*/false);
>>> +	// check if damping for this body is activated
>>> +	Real dampcoeff=(state->isDamped ? damping : 0);
>>> +	// always positive dissipation, by-component: |F_i|*|v_i|*dampcoeff*dt
>>> (| T_i|*|ω_i|*dampcoeff*dt for rotations)
>>> +	if(dampcoeff!=0.){
>>> +		scene->energy-
>>>
>>>> add(fluctVel.cwise().abs().dot(f.cwise().abs())*dampcoeff*scene-
>>>> dt,"nonviscDamp",nonviscDampIx,/*non-incremental*/false);
>>>>
>>>   		// when the aspherical integrator is used, torque is damped 
> instead
>>>   		of
>>>
>>> ang acceleration; this code is only approximate
>>> -		scene->energy->add(state-
>>>
>>>> angVel.cwise().abs().dot(m.cwise().abs())*damping*scene-
>>>> dt,"nonviscDamp",nonviscDampIx,false);
>>> +		scene->energy->add(state-
>>>
>>>> angVel.cwise().abs().dot(m.cwise().abs())*dampcoeff*scene-
>>>> dt,"nonviscDamp",nonviscDampIx,false);
>>>>
>>>   	}
>>>   	// kinetic energy
>>>   	Real Etrans=.5*state->mass*fluctVel.squaredNorm();
>>>
>>> @@ -106,9 +109,13 @@
>>>
>>>   	YADE_PARALLEL_FOREACH_BODY_BEGIN(const shared_ptr<Body>&  b,
>>>   	scene->bodies)
>>>
>>> {
>>>
>>>   			// clump members are handled inside clumps
>>>   			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();
>>>
>>> +
>>> +			// check if damping for this body is activated
>>> +			Real dampcoeff=(state->isDamped ? damping : 0);
>>> +
>>>
>>>   			// clumps forces
>>>   			if(b->isClump()) {
>>>   			
>>>   				b->shape-
>>>> cast<Clump>().addForceTorqueFromMembers(state,scene,f,m);
>>> @@ -146,7 +153,7 @@
>>>
>>>   			if (state->blockedDOFs!=State::DOF_ALL) {
>>>   			
>>>   				// linear acceleration
>>>   				Vector3r linAccel=computeAccel(f,state->mass,state-
>> blockedDOFs);
>>> -				cundallDamp2nd(dt,f,fluctVel,linAccel);
>>> +				cundallDamp2nd(dt,f,fluctVel,linAccel,dampcoeff);
>>>
>>>   				//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
>>>
>>> @@ -154,11 +161,11 @@
>>>
>>>   				// angular acceleration
>>>   				if(!useAspherical){ // uses angular velocity
>>>   				
>>>   					Vector3r angAccel=computeAngAccel(m,state->inertia,state-
>>>> blockedDOFs);
>>> -					cundallDamp2nd(dt,m,state->angVel,angAccel);
>>> +					cundallDamp2nd(dt,m,state->angVel,angAccel,dampcoeff);
>>>
>>>   					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
>>> -					cundallDamp1st(m,state->angVel);
>>> +					cundallDamp1st(m,state->angVel,dampcoeff);
>>>
>>>   				}
>>>   			
>>>   			}
>>>
>>> === modified file pkg/dem/NewtonIntegrator.hpp
>>> --- pkg/dem/NewtonIntegrator.hpp	2011-09-20 10:58:18 +0000
>>> +++ pkg/dem/NewtonIntegrator.hpp	2012-01-12 05:26:05 +0000
>>> @@ -22,8 +22,8 @@
>>>
>>>   class VelocityBins;
>>>   
>>>   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 cundallDamp1st(Vector3r&  force, const Vector3r&  vel,
>>> const Real&  dampcoeff);
>>> +	inline void cundallDamp2nd(const Real&  dt, const Vector3r&  force,
>>> const Vector3r&  vel, Vector3r&  accel, const Real&  dampcoeff);
>>>
>>>   	bool haveBins;
>>>   	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
>>>
>>> On Fri, 20 Jan 2012 02:06:51 AM Bruno Chareyre wrote:
>>>> Klaus,
>>>> If what you want is to increase the timestep (because your steel-net
>>>> induces very small dt, right?), this change will not work. You need to
>>>> modify density in the good place in Newton, not damping. It can be done
>>>> by adding dscale parameter in the material class (hence no need to check
>>>> if it exists), which would be 1 by default, hence no need to check if it
>>>> is NaN.
>>>> isnan is not in C++ standard, by the way. I know it is used in many
>>>> places but we should avoid it if possible. I'm also not sure of the cost
>>>> of this function, it may be more expensive than comparing two doubles.
>>>>
>>>> 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
>>> _______________________________________________
>>> 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
>> _______________________________________________
>> 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
> _______________________________________________
> 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
11, rue des Mathématiques
BP 46
38402 St Martin d'Hères, France
Tél : +33 4 56 52 86 21
Fax : +33 4 76 82 70 43
________________