yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #08228
Re: : Add damping coefficient to material
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.thoeni@xxxxxxxxx> 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.thoeni@xxxxxxxxx>
> 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@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
Follow ups
References
