One note: if(!isnan(b->getDampCoeff())) dampcoeff=b->getDampCoeff();
I think it is not optimal to check (if(!isnan(b->getDampCoeff()))) on _every_ interaction on _each_ step (performance!!!). As I understand, you do not use "damping"-variable of Newtonintegrator at all, so just set it to -1 and use it as flag. So you will check it 1 time per step. But it needs to be documented. I dont see any large problems with this commit. If you really need it, you can do it after fixing this issue. Anton On Thu, Jan 19, 2012 at 1:48 PM, Klaus Thoeni <klaus.tho...@gmail.com> wrote: > Here is the diff, effected classes Body, Material, NewtonIntegrator. Let me > know > if you are happy! > > === modified file core/Body.hpp > --- core/Body.hpp 2010-12-22 10:55:23 +0000 > +++ core/Body.hpp 2012-01-19 11:40:50 +0000 > @@ -68,6 +68,8 @@ > > int getGroupMask() {return groupMask; }; > bool maskOk(int mask){return (mask==0 || (groupMask&mask));} > + > + Real getDampCoeff() { return (!&Body::isClump) ? > material->damping : NaN; > }; > > // only BodyContainer can set the id of a body > friend class BodyContainer; > > === modified file core/Material.hpp > --- core/Material.hpp 2010-11-07 11:46:20 +0000 > +++ core/Material.hpp 2012-01-12 05:50:12 +0000 > @@ -39,7 +39,8 @@ > YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Material,Serializable,"Material > properties of a :yref:`body<Body>`.", > ((int,id,((void)"not shared",-1),Attr::readonly,"Numeric id of > this > material; is non-negative only if this Material is shared (i.e. in > O.materials), -1 otherwise. This value is set automatically when the material > is inserted to the simulation via > :yref:`O.materials.append<MaterialContainer.append>`. (This id was necessary > since before boost::serialization was used, shared pointers were not tracked > properly; it might disappear in the future)")) > ((string,label,,,"Textual identifier for this material; can be > used for > shared materials lookup in :yref:`MaterialContainer`.")) > - ((Real,density,1000,,"Density of the material [kg/m³]")), > + ((Real,density,1000,,"Density of the material [kg/m³]")) > + ((Real,damping,NaN,,"Local damping coefficient associated > with the > material [-]. If this value is given :yref:`Newtonintegrator` uses this value > instead of :yref:`damping<NewtonIntegrator.damping>` and it allows to use > different damping coefficients for different materials.")), > /* ctor */, > /*py*/ > .def("newAssocState",&Material::newAssocState,"Return new > :yref:`State` > instance, which is associated with this :yref:`Material`. Some materials have > special requirement on :yref:`Body::state` type and calling this function when > the body is created will ensure that they match. (This is done automatically > if you use utils.sphere, … functions from python).") > > === modified file pkg/dem/NewtonIntegrator.cpp > --- pkg/dem/NewtonIntegrator.cpp 2011-11-30 17:39:33 +0000 > +++ pkg/dem/NewtonIntegrator.cpp 2012-01-19 11:56:41 +0000 > @@ -11,17 +11,18 @@ > #include<yade/pkg/dem/Clump.hpp> > #include<yade/pkg/common/VelocityBins.hpp> > #include<yade/lib/base/Math.hpp> > +#include "../../lib/base/Logging.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,14 @@ > > 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 which damping coefficient to use > + Real dampcoeff=damping; > + if(!isnan(b->getDampCoeff())) dampcoeff=b->getDampCoeff(); > + // 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(); > @@ -107,6 +111,10 @@ > // clump members are handled inside clumps > if(unlikely(b->isClumpMember())) continue; > > + // check which damping coefficient to use > + Real dampcoeff=damping; > + if(!isnan(b->getDampCoeff())) > dampcoeff=b->getDampCoeff(); > + > State* state=b->state.get(); const Body::id_t& > id=b->getId(); > Vector3r f=gravity*state->mass, m=Vector3r::Zero(); > // clumps forces > @@ -146,7 +154,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 +162,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 Thu, 19 Jan 2012 11:18:27 PM Anton Gladky wrote: >> Could you, please, attach a diff first? >> Thanks >> >> Anton >> >> On Thu, Jan 19, 2012 at 1:03 PM, Klaus Thoeni <klaus.tho...@gmail.com> > wrote: >> > Solved the problem by checking if it is a clump: >> > >> > Real getDampCoeff() { return (!&Body::isClump) ? material->damping : NaN; >> > }; >> > >> > However, I think we still have to sort out why clumps don't have a >> > material! >> > >> > I will commit the code, it might be useful to have different non-viscous >> > damping coefficients, or what do you think? >> > >> > Klaus >> >> _______________________________________________ >> 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