yade-dev team mailing list archive

[Klaus Thoeni <klaus.thoeni@xxxxxxxxx>]

Hi guys,

thanks a lot for your suggestions. I just committed the code and I think it 
looks very good now.

Klaus

On Mon, 23 Jan 2012 11:08:59 PM Bruno Chareyre wrote:
> 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