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