yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #08289
[Branch ~yade-dev/yade/trunk] Rev 3017: - Compute the critical timestep correctly when clumps are present.
------------------------------------------------------------
revno: 3017
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
branch nick: trunk
timestamp: Thu 2012-02-02 18:56:00 +0100
message:
- Compute the critical timestep correctly when clumps are present.
modified:
pkg/dem/GlobalStiffnessTimeStepper.cpp
pkg/dem/GlobalStiffnessTimeStepper.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/GlobalStiffnessTimeStepper.cpp'
--- pkg/dem/GlobalStiffnessTimeStepper.cpp 2011-11-24 14:20:42 +0000
+++ pkg/dem/GlobalStiffnessTimeStepper.cpp 2012-02-02 17:56:00 +0000
@@ -12,6 +12,7 @@
#include<yade/pkg/dem/DemXDofGeom.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/Scene.hpp>
+#include<yade/pkg/dem/Clump.hpp>
CREATE_LOGGER(GlobalStiffnessTimeStepper);
YADE_PLUGIN((GlobalStiffnessTimeStepper));
@@ -20,16 +21,23 @@
void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& body, Scene * ncb)
{
- const State* sdec=body->state.get();
-
+ const 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(!sdec || stiffness==Vector3r::Zero())
return; // not possible to compute!
- Real dtx, dty, dtz;
-
+ 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
@@ -37,12 +45,10 @@
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;
-
+ 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;
@@ -69,21 +75,16 @@
// shared_ptr<InteractionContainer>& interactions = ncb->interactions;
newDt = Mathr::MAX_REAL;
+ computedSomething=false;
BodyContainer::iterator bi = bodies->begin();
BodyContainer::iterator biEnd = bodies->end();
- for( ; bi!=biEnd ; ++bi )
- {
-// if (!*bi) continue;
+ for( ; bi!=biEnd ; ++bi ){
shared_ptr<Body> b = *bi;
- if (b->isDynamic()) findTimeStepFromBody(b, ncb);
- }
-
- if(computedSomething)
- {
- previousDt = min ( min(newDt , defaultDt), 1.5*previousDt );// at maximum, dt will be multiplied by 1.5 in one iterration, this is to prevent brutal switches from 0.000... to 1 in some computations
+ if (b->isDynamic() && !b->isClumpMember()) findTimeStepFromBody(b, ncb);}
+ if(computedSomething){
+ 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;
- }
+ computedOnce = true;}
else if (!computedOnce) scene->dt=defaultDt;
LOG_INFO("computed timestep " << newDt <<
(scene->dt==newDt ? string(", appplied") :
=== modified file 'pkg/dem/GlobalStiffnessTimeStepper.hpp'
--- pkg/dem/GlobalStiffnessTimeStepper.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/GlobalStiffnessTimeStepper.hpp 2012-02-02 17:56:00 +0000
@@ -38,7 +38,8 @@
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",
- ((Real,defaultDt,1,,"used as default AND as max value of the timestep"))
+ ((Real,defaultDt,1,,"used as the initial value of the timestep (especially usefull in the first steps when no contact exist). Can be defined using :yref:`utils.PWaveTimeStep`"))
+ ((Real,maxDt,1,,"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))")),
computedOnce=false;)
