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@lists.launchpad.net
> >> Unsubscribe : https://launchpad.net/~yade-dev
> >> More help : https://help.launchpad.net/ListHelp
> >
> > _______________________________________________
> > Mailing list: https://launchpad.net/~yade-dev
> > Post to : yade-dev@lists.launchpad.net
> > Unsubscribe : https://launchpad.net/~yade-dev
> > More help : https://help.launchpad.net/ListHelp
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to : yade-dev@lists.launchpad.net
> Unsubscribe : https://launchpad.net/~yade-dev
> More help : https://help.launchpad.net/ListHelp
_______________________________________________
Mailing list: https://launchpad.net/~yade-dev
Post to : yade-dev@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-dev
More help : https://help.launchpad.net/ListHelp
