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 <[email protected]>: > 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 > À: [email protected] > 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 * > * [email protected] * > * * > * 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 * > * [email protected] * > * * > * 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 : [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
