yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #09995
Re: Modified GlobalStiffnessTimeStepper for visco-elastic contact law
Hi Raphaël,
could you, please, attach a diff (or patch) instead of copying it into
the mail-body.
Thanks,
Anton
2013/9/20 Raphaël Maurin <raphael.maurin@xxxxxxxxx>:
> Hi all,
> Here is the final version of GlobalStiffnessTimeStepper that I modified to
> take into account the viscous contribution when using
> Law2_ScGeom_ViscElPhys_Basic.
> I didn't modify the elastic part, so that it will work exactly the same for
> people using it as before (considering only elastic contribution). I just
> added an option : viscEl which, when True, allows to calculate the time step
> for ViscElPhys as described before.
> I tested it and everything works well.
>
> I put the patch below for GlobalStiffnessTimeStepper.cpp and
> GlobalStiffnessTimeStepper.hpp, can someone add the change on github ?
> Thanks
>
> Raphael
>
>
>
> GlobalStiffnessTimeStepper.cpp
>
>
>
> diff --git a/pkg/dem/GlobalStiffnessTimeStepper.cpp
> b/pkg/dem/GlobalStiffnessTimeStepper.cpp
> index 8a41ceb..0205655 100644
> --- a/pkg/dem/GlobalStiffnessTimeStepper.cpp
> +++ b/pkg/dem/GlobalStiffnessTimeStepper.cpp
> @@ -14,6 +14,7 @@
> #include<yade/core/Scene.hpp>
> #include<yade/core/Clump.hpp>
> #include<yade/pkg/dem/Shop.hpp>
> +#include<yade/pkg/dem/ViscoelasticPM.hpp>
>
> CREATE_LOGGER(GlobalStiffnessTimeStepper);
> YADE_PLUGIN((GlobalStiffnessTimeStepper));
> @@ -25,13 +26,16 @@ void
> GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
> State* sdec=body->state.get();
> Vector3r& stiffness= stiffnesses[body->getId()];
> Vector3r& Rstiffness=Rstiffnesses[body->getId()];
> -
> if(body->isClump()) {// if clump, we sum stifnesses of all members
> const shared_ptr<Clump>&
> clump=YADE_PTR_CAST<Clump>(body->shape);
> FOREACH(Clump::MemberMap::value_type& B, clump->members){
> const shared_ptr<Body>& b =
> Body::byId(B.first,scene);
> stiffness+=stiffnesses[b->getId()];
> Rstiffness+=Rstiffnesses[b->getId()];
> + if (viscEl == true){
> +
> viscosities[body->getId()]+=viscosities[b->getId()];
> +
> Rviscosities[body->getId()]+=Rviscosities[b->getId()];
> + }
> }
> }
>
> @@ -39,13 +43,15 @@ void
> GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
> if (densityScaling) sdec->densityScaling =
> min(1.0001*sdec->densityScaling,
> timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
> return; // not possible to compute!
> }
> -
> +
> + //Determine the elastic minimum eigenperiod (and if required
> determine also the viscous one separately and take the minimum of the two)
> +
> + //Elastic
> Real dtx, dty, dtz;
> Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
> if (dt!=0) {
> dt = sdec->mass/dt; computedSomething = true;}//dt =
> squared eigenperiod of translational motion
> else dt = Mathr::MAX_REAL;
> -
> if (Rstiffness.x()!=0) {
> dtx = sdec->inertia.x()/Rstiffness.x(); computedSomething =
> true;}//dtx = squared eigenperiod of rotational motion around x
> else dtx = Mathr::MAX_REAL;
> @@ -56,15 +62,43 @@ void
> GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& bo
> dtz = sdec->inertia.z()/Rstiffness.z(); computedSomething =
> true;}
> else dtz = Mathr::MAX_REAL;
>
> - Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest
> squared eigenperiod for rotational motions
> + Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest
> squared eigenperiod for elastic rotational motions
> dt =
> 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 =
> sqrt(2)
> +
> + //Viscous
> + if (viscEl == true){
> + Vector3r& viscosity = viscosities[body->getId()];
> + Vector3r& Rviscosity = Rviscosities[body->getId()];
> + Real dtx_visc, dty_visc, dtz_visc;
> + Real dt_visc = max(max(viscosity.x(), viscosity.y()),
> viscosity.z() );
> + if (dt_visc!=0) {
> + dt_visc = sdec->mass/dt_visc; computedSomething =
> true;}//dt = eigenperiod of the viscous translational motion
> + else {dt_visc = Mathr::MAX_REAL;}
> +
> + if (Rviscosity.x()!=0) {
> + dtx_visc = sdec->inertia.x()/Rviscosity.x();
> computedSomething = true;}//dtx = eigenperiod of viscous rotational motion
> around x
> + else dtx_visc = Mathr::MAX_REAL;
> + if (Rviscosity.y()!=0) {
> + dty_visc = sdec->inertia.y()/Rviscosity.y();
> computedSomething = true;}
> + else dty_visc = Mathr::MAX_REAL;
> + if (Rviscosity.z()!=0) {
> + dtz_visc = sdec->inertia.z()/Rviscosity.z();
> computedSomething = true;}
> + else dtz_visc = Mathr::MAX_REAL;
> +
> + Real Rdt_visc = std::min( std::min (dtx_visc, dty_visc),
> dtz_visc );//Rdt = smallest squared eigenperiod for viscous rotational
> motions
> + dt_visc =
> 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);
> +
> + //Take the minimum between the elastic and viscous minimum
> eigenperiod.
> + dt = std::min(dt,dt_visc);
> + }
> +
> //if there is a target dt, then we apply density scaling on the
> body, the inertia used in Newton will be mass*scaling, the weight is
> unmodified
> if (densityScaling) {
> sdec->densityScaling =
> min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
> newDt=targetDt;
> }
> //else we update dt normaly
> - else {newDt = std::min(dt,newDt);}
> + else {newDt = std::min(dt,newDt);}
> }
>
> bool GlobalStiffnessTimeStepper::isActivated()
> @@ -97,7 +131,7 @@ void GlobalStiffnessTimeStepper::computeTimeStep(Scene*
> ncb)
> computedOnce = true;}
> else if (!computedOnce) scene->dt=defaultDt;
> LOG_INFO("computed timestep " << newDt <<
> - (scene->dt==newDt ? string(", appplied") :
> + (scene->dt==newDt ? string(", applied") :
> string(", BUT timestep is
> ")+lexical_cast<string>(scene->dt))<<".");
> }
>
> @@ -107,15 +141,24 @@ void
> GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
> if(size<rb->bodies->size()){
> size=rb->bodies->size();
> stiffnesses.resize(size); Rstiffnesses.resize(size);
> + if (viscEl == true){
> + viscosities.resize(size); Rviscosities.resize(size);
> + }
> }
> /* reset stored values */
> memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
> memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
> + if (viscEl == true){
> + memset(& viscosities[0],0,sizeof(Vector3r)*size);
> + memset(&Rviscosities[0],0,sizeof(Vector3r)*size);
> + }
> +
> FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
> if(!contact->isReal()) continue;
>
> GenericSpheresContact*
> geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
> NormShearPhys*
> phys=YADE_CAST<NormShearPhys*>(contact->phys.get()); assert(phys);
> +
> // all we need for getting stiffness
> Vector3r& normal=geom->normal; Real& kn=phys->kn; Real&
> ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
> Real fn = (static_cast<NormShearPhys *>
> (contact->phys.get()))->normalForce.squaredNorm();
> @@ -134,10 +177,35 @@ void
> GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
>
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> diag_Rstiffness *= ks;
> +
> +
> //NOTE : contact laws with moments would be handled
> correctly by summing directly bending+twisting stiffness to diag_Rstiffness.
> The fact that there is no problem currently with e.g. cohesiveFrict law is
> probably because final computed dt is constrained by translational motion,
> not rotations.
> stiffnesses [contact->getId1()]+=diag_stiffness;
>
> Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
> stiffnesses [contact->getId2()]+=diag_stiffness;
> -
> Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
> +
> Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
> +
> + //Same for the Viscous part, if required
> + if (viscEl == true){
> + ViscElPhys* viscPhys =
> YADE_CAST<ViscElPhys*>(contact->phys.get()); assert(viscPhys);
> + Real& cn=viscPhys->cn; Real& cs=viscPhys->cs;
> + //Diagonal terms of the translational viscous matrix
> + Vector3r diag_viscosity =
> Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
> + diag_viscosity *= cn-cs;
> + diag_viscosity = diag_viscosity +
> Vector3r(1,1,1)*cs;
> + //diagonal terms of the rotational viscous matrix
> + Vector3r diag_Rviscosity =
> +
> Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
> +
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
> +
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> + diag_Rviscosity *= cs;
> +
> + // Add the contact stiffness matrix to the two
> particles one
> + viscosities [contact->getId1()]+=diag_viscosity;
> +
> Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
> + viscosities [contact->getId2()]+=diag_viscosity;
> +
> Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
> + }
> +
> }
> }
>
>
>
>
>
>
> GlobalStiffnessTimeStepper.hpp
>
>
>
>
>
>
> diff --git a/pkg/dem/GlobalStiffnessTimeStepper.hpp
> b/pkg/dem/GlobalStiffnessTimeStepper.hpp
> index 1d60f88..bd657fa 100644
> --- a/pkg/dem/GlobalStiffnessTimeStepper.hpp
> +++ b/pkg/dem/GlobalStiffnessTimeStepper.hpp
> @@ -24,6 +24,8 @@ class GlobalStiffnessTimeStepper : public TimeStepper
> private :
> vector<Vector3r> stiffnesses;
> vector<Vector3r> Rstiffnesses;
> + vector<Vector3r> viscosities;
> + vector<Vector3r> Rviscosities;
> void computeStiffnesses(Scene*);
>
> Real newDt;
> @@ -37,13 +39,14 @@ class GlobalStiffnessTimeStepper : public TimeStepper
> virtual void computeTimeStep(Scene*);
> virtual bool isActivated();
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(
> - GlobalStiffnessTimeStepper,TimeStepper,"An engine
> assigning the time-step as a fraction of the minimum eigen-period in the
> problem. The derivation is detailed in the chapter on `DEM formulation
> <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
> + GlobalStiffnessTimeStepper,TimeStepper,"An engine
> assigning the time-step as a fraction of the minimum eigen-period in the
> problem. The derivation is detailed in the chapter on `DEM formulation
> <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_. The
> viscEl option enables to evaluate the timestep in a similar way for the
> visco-elastic contact law :yref:`Law2_ScGeom_ViscElPhys_Basic`, more detail
> in :yref:`GlobalStiffnessTimestepper::viscEl`. ",
> ((Real,defaultDt,-1,,"used as the initial value of
> the timestep (especially useful in the first steps when no contact exist).
> If negative, it will be defined by :yref:`utils.PWaveTimeStep` *
> :yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
> ((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as
> max value of the timestep whatever the computed value"))
> ((Real,previousDt,1,,"last computed dt |yupdate|"))
> ((Real,timestepSafetyCoefficient,0.8,,"safety factor
> between the minimum eigen-period and the final assigned dt (less than 1))"))
> ((bool,densityScaling,false,,"|yupdate| don't modify
> this value if you don't plan to modify the scaling factor manually for some
> bodies. In most cases, it is enough to set
> :yref:`NewtonIntegrator::densityScaling` and let this one be adjusted
> automatically."))
> - ((Real,targetDt,1,,"if
> :yref:`NewtonIntegrator::densityScaling` is active, this value will be used
> as the simulation timestep and the scaling will use this value of dt as the
> target value. The value of targetDt is arbitrary and should have no effect
> in the result in general. However if some bodies have imposed velocities,
> for instance, they will move more or less per each step depending on this
> value.")),
> + ((Real,targetDt,1,,"if
> :yref:`NewtonIntegrator::densityScaling` is active, this value will be used
> as the simulation timestep and the scaling will use this value of dt as the
> target value. The value of targetDt is arbitrary and should have no effect
> in the result in general. However if some bodies have imposed velocities,
> for instance, they will move more or less per each step depending on this
> value."))
> + ((bool,viscEl,false,,"To use with
> :yref:`ViscElPhys`. if True, evaluate separetly the minimum eigen-period in
> the problem considering only the elastic contribution on one hand (spring
> only), and only the viscous contribution on the other hand (dashpot only).
> Take then the minimum of the two and use the safety coefficient
> :yref:`GlobalStiffnessTimestepper::timestepSafetyCoefficient` to take into
> account the possible coupling between the two contribution.")),
> computedOnce=false;)
> DECLARE_LOGGER;
> };
>
>
>
>
>
>
>
>
>
>
> -------- Message d'origine--------
> De: Raphaël Maurin
> Date: mar. 03/09/2013 12:21
> À: yade-dev@xxxxxxxxxxxxxxxxxxx
> Objet : Modified GlobalStiffnessTimeStepper for visco-elastic contact law
>
> Hi all,
>
> I developed a "copy" of GlobalStiffnessTimeStepper which evaluate the time
> step of a system of particle with Law2_ScGeom_ViscElPhys_Basic contact law
> (spring + viscous damping).
>
> The principle of the time step evaluation is :
> We consider that the elastic and the viscous degrees of freedom are
> decoupled. From that, we evaluate the critical time step of the system of
> particles considering the elastic part only (dt_el) on one hand, and
> considering the viscous part only (dt_visc) on the other hand. Then, we take
> the minimum of the two timestep and multiply it by a safety factor to
> account for the possible coupling between the two effect (dt =
> alpha*min(dt_el,dt_visc).
>
> To do so, I took GlobalStiffnessTimeStepper (GSTS) which is already
> evaluating the elastic critical time step, and implemented similarly a
> calculation for the viscous critical time step of the system from the the
> viscosity matrix of each particles. The changes are not so important between
> the version I implemented and the last one (see the code below). I tested it
> in cases where the normal viscous part is dominant and the new version gives
> a lower time step which allows to keep the caculus stable, while in the old
> version everything explode, so it seems to work.
>
> I have some questions about that :
> - Does the evaluation of the time step seems reasonable to you?
> - What can I do to test it quantitavely ? (if I try a collision between two
> beads, GSTS does not work..)
> - To add it to the code, should I add it as a new "file"/"function" or
> should I implement it as an option of GSTS ?
>
> And about GSTS in general :
> - why is GSTS not working when we consider a gravity deposition or a
> collision ?? (even when defaultDt is small enough)
>
> Thanks,
>
> Raphael
>
>
>
> GlobalStiffnessTimeStepper.cpp
>
> /*************************************************************************
> * Copyright (C) 2006 by Bruno Chareyre *
> * bruno.chareyre@xxxxxxxxxxx *
> * *
> * This program is free software; it is licensed under the terms of the *
> * GNU General Public License v2 or later. See file LICENSE for details. *
> *************************************************************************/
>
> #include"GlobalStiffnessTimeStepper.hpp"
> #include<yade/pkg/dem/FrictPhys.hpp>
> #include<yade/pkg/dem/ScGeom.hpp>
> #include<yade/pkg/dem/DemXDofGeom.hpp>
> #include<yade/core/Interaction.hpp>
> #include<yade/core/Scene.hpp>
> #include<yade/core/Clump.hpp>
> #include<yade/pkg/dem/Shop.hpp>
> #include<yade/pkg/dem/ViscoelasticPM.hpp>
>
> CREATE_LOGGER(GlobalStiffnessTimeStepper);
> YADE_PLUGIN((GlobalStiffnessTimeStepper));
>
> GlobalStiffnessTimeStepper::~GlobalStiffnessTimeStepper() {}
>
> void GlobalStiffnessTimeStepper::findTimeStepFromBody(const
> shared_ptr<Body>& body, Scene * ncb)
> {
> State* sdec=body->state.get();
> Vector3r& stiffness= stiffnesses[body->getId()];
> Vector3r& Rstiffness=Rstiffnesses[body->getId()];
> Vector3r& viscosity= viscosities[body->getId()];
> Vector3r& Rviscosity= Rviscosities[body->getId()];
>
> if(body->isClump()) {// if clump, we sum stifnesses of all members
> const shared_ptr<Clump>&
> clump=YADE_PTR_CAST<Clump>(body->shape);
> FOREACH(Clump::MemberMap::value_type& B, clump->members){
> const shared_ptr<Body>& b =
> Body::byId(B.first,scene);
> stiffness+=stiffnesses[b->getId()];
> Rstiffness+=Rstiffnesses[b->getId()];
> viscosity+=viscosities[b->getId()];
> Rviscosity+=Rviscosities[b->getId()];
> }
> }
>
> if(!sdec || stiffness==Vector3r::Zero()){
> if (densityScaling) sdec->densityScaling =
> min(1.0001*sdec->densityScaling,
> timestepSafetyCoefficient*pow(defaultDt/targetDt,2.0));
> return; // not possible to compute!
> }
>
> //Determine the elastic critical time step and the viscous one
> separately (as if they were decoupled) and take the minimum time step
> with a safety coefficient.
>
> //Elastic
> Real dtx, dty, dtz;
> Real dt = max( max (stiffness.x(), stiffness.y()), stiffness.z() );
> if (dt!=0) {
> dt = sdec->mass/dt; computedSomething = true;}//dt =
> squared eigenperiod of translational motion
> else dt = Mathr::MAX_REAL;
> if (Rstiffness.x()!=0) {
> dtx = sdec->inertia.x()/Rstiffness.x(); computedSomething =
> true;}//dtx = squared eigenperiod of rotational motion around x
> else dtx = Mathr::MAX_REAL;
> if (Rstiffness.y()!=0) {
> dty = sdec->inertia.y()/Rstiffness.y(); computedSomething =
> true;}
> else dty = Mathr::MAX_REAL;
> if (Rstiffness.z()!=0) {
> dtz = sdec->inertia.z()/Rstiffness.z(); computedSomething =
> true;}
> else dtz = Mathr::MAX_REAL;
>
> Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest
> squared eigenperiod for elastic rotational motions
> dt =
> 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 =
> sqrt(2)
>
> //Viscous
> Real dtx_visc, dty_visc, dtz_visc;
> Real dt_visc = max(max(viscosity.x(), viscosity.y()), viscosity.z()
> );
> if (dt_visc!=0) {
> dt_visc = sdec->mass/dt_visc; computedSomething =
> true;}//dt = eigenperiod of the viscous translational motion
> else {dt_visc = Mathr::MAX_REAL;}
>
> if (Rviscosity.x()!=0) {
> dtx_visc = sdec->inertia.x()/Rviscosity.x();
> computedSomething = true;}//dtx = eigenperiod of viscous rotational motion
> around x
> else dtx_visc = Mathr::MAX_REAL;
> if (Rviscosity.y()!=0) {
> dty_visc = sdec->inertia.y()/Rviscosity.y();
> computedSomething = true;}
> else dty_visc = Mathr::MAX_REAL;
> if (Rviscosity.z()!=0) {
> dtz_visc = sdec->inertia.z()/Rviscosity.z();
> computedSomething = true;}
> else dtz_visc = Mathr::MAX_REAL;
>
> Real Rdt_visc = std::min( std::min (dtx_visc, dty_visc), dtz_visc
> );//Rdt = smallest squared eigenperiod for viscous rotational motions
> dt_visc = 2*timestepSafetyCoefficient*std::min(dt_visc,Rdt_visc);
>
> //Take the minimum between the elastic and viscous critical time
> step.
> dt = std::min(dt,dt_visc);
>
> //if there is a target dt, then we apply density scaling on the
> body, the inertia used in Newton will be mass*scaling, the weight is
> unmodified
> if (densityScaling) {
> sdec->densityScaling =
> min(sdec->densityScaling,timestepSafetyCoefficient*pow(dt /targetDt,2.0));
> newDt=targetDt;
> }
> //else we update dt normaly
> else {newDt = std::min(dt,newDt);}
> }
>
> bool GlobalStiffnessTimeStepper::isActivated()
> {
> return (active && ((!computedOnce) || (scene->iter %
> timeStepUpdateInterval == 0) || (scene->iter < (long int) 2) ));
> }
>
> void GlobalStiffnessTimeStepper::computeTimeStep(Scene* ncb)
> {
> // for some reason, this line is necessary to have correct
> functioning (no idea _why_)
> // see scripts/test/compare-identical.py, run with or without
> active=active.
> active=active;
> if (defaultDt<0) defaultDt=
> timestepSafetyCoefficient*Shop::PWaveTimeStep(Omega::instance().getScene());
> computeStiffnesses(ncb);
>
> shared_ptr<BodyContainer>& bodies = ncb->bodies;
> newDt = Mathr::MAX_REAL;
> computedSomething=false;
> BodyContainer::iterator bi = bodies->begin();
> BodyContainer::iterator biEnd = bodies->end();
> for( ; bi!=biEnd ; ++bi ){
> shared_ptr<Body> b = *bi;
> if (b->isDynamic() && !b->isClumpMember())
> findTimeStepFromBody(b, ncb);
>
> }
> if(densityScaling) (newDt=targetDt);
> if(computedSomething || densityScaling){
> previousDt = min ( min(newDt , maxDt), 1.05*previousDt );//
> at maximum, dt will be multiplied by 1.05 in one iterration, this is to
> prevent brutal switches from 0.000... to 1 in some computations
> scene->dt=previousDt;
> computedOnce = true;}
> else if (!computedOnce) scene->dt=defaultDt;
> LOG_INFO("computed timestep " << newDt <<
> (scene->dt==newDt ? string(", applied") :
> string(", BUT timestep is
> ")+lexical_cast<string>(scene->dt))<<".");
> }
>
> void GlobalStiffnessTimeStepper::computeStiffnesses(Scene* rb){
> /* check size */
> size_t size=stiffnesses.size();
> if(size<rb->bodies->size()){
> size=rb->bodies->size();
> stiffnesses.resize(size); Rstiffnesses.resize(size);
> viscosities.resize(size); Rviscosities.resize(size);
> }
> /* reset stored values */
> memset(& stiffnesses[0],0,sizeof(Vector3r)*size);
> memset(&Rstiffnesses[0],0,sizeof(Vector3r)*size);
> memset(& viscosities[0],0,sizeof(Vector3r)*size);
> memset(& Rviscosities[0],0,sizeof(Vector3r)*size);
>
> FOREACH(const shared_ptr<Interaction>& contact, *rb->interactions){
> if(!contact->isReal()) continue;
>
> GenericSpheresContact*
> geom=YADE_CAST<GenericSpheresContact*>(contact->geom.get()); assert(geom);
> ViscElPhys* phys =
> YADE_CAST<ViscElPhys*>(contact->phys.get()); assert(phys);
> // all we need for getting stiffness
> Vector3r& normal=geom->normal; Real& kn=phys->kn; Real&
> ks=phys->ks; Real& radius1=geom->refR1; Real& radius2=geom->refR2;
> Real& cn=phys->cn; Real& cs=phys->cs;
> Real fn = (static_cast<NormShearPhys *>
> (contact->phys.get()))->normalForce.squaredNorm();
> if (fn==0) continue;//Is it a problem with some laws? I
> still don't see why.
>
> //Diagonal terms of the translational stiffness matrix
> Vector3r diag_stiffness =
> Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
> diag_stiffness *= kn-ks;
> diag_stiffness = diag_stiffness + Vector3r(1,1,1)*ks;
>
> //diagonal terms of the rotational stiffness matrix
> // Vector3r branch1 =
> currentContactGeometry->normal*currentContactGeometry->radius1;
> // Vector3r branch2 =
> currentContactGeometry->normal*currentContactGeometry->radius2;
> Vector3r diag_Rstiffness =
>
> Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> diag_Rstiffness *= ks;
>
> //Diagonal terms of the translational viscous matrix
> Vector3r diag_viscosity =
> Vector3r(std::pow(normal.x(),2),std::pow(normal.y(),2),std::pow(normal.z(),2));
> diag_viscosity *= cn-cs;
> diag_viscosity = diag_viscosity + Vector3r(1,1,1)*cs;
> //diagonal terms of the rotational viscous matrix
> Vector3r diag_Rviscosity =
>
> Vector3r(std::pow(normal.y(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.z(),2),
>
> std::pow(normal.x(),2)+std::pow(normal.y(),2));
> diag_Rviscosity *= cs;
>
> //NOTE : contact laws with moments would be handled
> correctly by summing directly bending+twisting stiffness to diag_Rstiffness.
> The fact that there is no problem currently with e.g. cohesiveFrict law is
> probably because final computed dt is constrained by translational motion,
> not rotations.
> stiffnesses [contact->getId1()]+=diag_stiffness;
>
> Rstiffnesses[contact->getId1()]+=diag_Rstiffness*pow(radius1,2);
> viscosities [contact->getId1()]+=diag_viscosity;
>
> Rviscosities[contact->getId1()]+=diag_Rviscosity*pow(radius1,2);
> stiffnesses [contact->getId2()]+=diag_stiffness;
>
> Rstiffnesses[contact->getId2()]+=diag_Rstiffness*pow(radius2,2);
> viscosities [contact->getId2()]+=diag_viscosity;
>
> Rviscosities[contact->getId2()]+=diag_Rviscosity*pow(radius2,2);
> }
> }
>
>
>
>
>
>
>
>
>
>
>
> GlobalStiffnessTimeStepper.hpp
>
>
>
> /*************************************************************************
> * Copyright (C) 2006 by Bruno Chareyre *
> * bruno.chareyre@xxxxxxxxxxx *
> * *
> * This program is free software; it is licensed under the terms of the *
> * GNU General Public License v2 or later. See file LICENSE for details. *
> *************************************************************************/
>
> #pragma once
>
> #include<yade/core/TimeStepper.hpp>
>
>
> /*! \brief Compute the critical timestep of the leap-frog scheme based on
> global stiffness of bodies
> See usage details in TriaxialTest documentation (TriaxialTest is
> also a good example of how to use this class)
> */
>
> class Interaction;
> class BodyContainer;
> class Scene;
>
> class GlobalStiffnessTimeStepper : public TimeStepper
> {
> private :
> vector<Vector3r> stiffnesses;
> vector<Vector3r> Rstiffnesses;
> vector<Vector3r> viscosities;
> vector<Vector3r> Rviscosities;
> void computeStiffnesses(Scene*);
>
> Real newDt;
> bool computedSomething,
> computedOnce;
> void findTimeStepFromBody(const shared_ptr<Body>& body,
> Scene * ncb);
>
> public :
> virtual ~GlobalStiffnessTimeStepper();
>
> virtual void computeTimeStep(Scene*);
> virtual bool isActivated();
> YADE_CLASS_BASE_DOC_ATTRS_CTOR(
> GlobalStiffnessTimeStepper,TimeStepper,"An engine
> assigning the time-step as a fraction of the minimum eigen-period in the
> problem. The derivation is detailed in the chapter on `DEM formulation
> <https://www.yade-dem.org/doc/formulation.html#dem-simulations>`_",
> ((Real,defaultDt,-1,,"used as the initial value of
> the timestep (especially useful in the first steps when no contact exist).
> If negative, it will be defined by :yref:`utils.PWaveTimeStep` *
> :yref:`GlobalStiffnessTimeStepper::timestepSafetyCoefficient`"))
> ((Real,maxDt,Mathr::MAX_REAL,,"if positive, used as
> max value of the timestep whatever the computed value"))
> ((Real,previousDt,1,,"last computed dt |yupdate|"))
> ((Real,timestepSafetyCoefficient,0.8,,"safety factor
> between the minimum eigen-period and the final assigned dt (less than 1))"))
> ((bool,densityScaling,false,,"|yupdate| don't modify
> this value if you don't plan to modify the scaling factor manually for some
> bodies. In most cases, it is enough to set
> :yref:`NewtonIntegrator::densityScaling` and let this one be adjusted
> automatically."))
> ((Real,targetDt,1,,"if
> :yref:`NewtonIntegrator::densityScaling` is active, this value will be used
> as the simulation timestep and the scaling will use this value of dt as the
> target value. The value of targetDt is arbitrary and should have no effect
> in the result in general. However if some bodies have imposed velocities,
> for instance, they will move more or less per each step depending on this
> value.")),
> computedOnce=false;)
> DECLARE_LOGGER;
> };
>
> REGISTER_SERIALIZABLE(GlobalStiffnessTimeStepper);
>
>
>
>
>
> _______________________________________________
> 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
