yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03445
[Branch ~yade-dev/yade/trunk] Rev 2038: - Fix the bug (r2025 probably) in GSTimeStepper (test of line 28).
------------------------------------------------------------
revno: 2038
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Thu 2010-02-18 23:47:05 +0100
message:
- Fix the bug (r2025 probably) in GSTimeStepper (test of line 28).
- Register timeSteppers.
- numpy : comment and new line at EOF
- TW : fix a warning.
- PTCollider : remove a "using namespace std" in the header.
modified:
core/TimeStepper.cpp
core/TimeStepper.hpp
lib/pyutil/numpy.hpp
pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.hpp
pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp
pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp
pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp
--
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 'core/TimeStepper.cpp'
--- core/TimeStepper.cpp 2009-12-09 17:11:51 +0000
+++ core/TimeStepper.cpp 2010-02-18 22:47:05 +0000
@@ -12,11 +12,11 @@
#include<yade/core/GlobalEngine.hpp>
#include<yade/core/Scene.hpp>
-TimeStepper::TimeStepper() : GlobalEngine()
-{
- active = true;
- timeStepUpdateInterval = 1;
-}
+// TimeStepper::TimeStepper() : GlobalEngine()
+// {
+// active = true;
+// timeStepUpdateInterval = 1;
+// }
bool TimeStepper::isActivated(Scene* mb)
=== modified file 'core/TimeStepper.hpp'
--- core/TimeStepper.hpp 2009-12-09 17:11:51 +0000
+++ core/TimeStepper.hpp 2010-02-18 22:47:05 +0000
@@ -18,17 +18,22 @@
class TimeStepper : public GlobalEngine
{
public :
- bool active;
- unsigned int timeStepUpdateInterval;
+// bool active;
+// unsigned int timeStepUpdateInterval;
- TimeStepper();
+// TimeStepper();
virtual void computeTimeStep(Scene* ) { throw; };
virtual bool isActivated(Scene*);
virtual void action(Scene* b) { computeTimeStep(b);} ;
void setActive(bool a, int nb=-1);
-
+
+ YADE_CLASS_BASE_DOC_ATTRS(
+ TimeStepper,GlobalEngine,"Engine defining time-step (fundamental class)",
+ ((bool,active,true,"is the engine active?"))
+ ((unsigned int,timeStepUpdateInterval,1,"dt update interval")));
+/*
REGISTER_ATTRIBUTES(GlobalEngine,(active)(timeStepUpdateInterval));
- REGISTER_CLASS_AND_BASE(TimeStepper,GlobalEngine);
+ REGISTER_CLASS_AND_BASE(TimeStepper,GlobalEngine);*/
};
REGISTER_SERIALIZABLE(TimeStepper);
=== modified file 'lib/pyutil/numpy.hpp'
--- lib/pyutil/numpy.hpp 2010-01-21 16:02:38 +0000
+++ lib/pyutil/numpy.hpp 2010-02-18 22:47:05 +0000
@@ -2,7 +2,6 @@
#pragma once
#include"numpy_boost.hpp"
-// helper macro do assign Vector3r values to a subarray
+// helper macro do assign Vector3r and Matrix3r values to subarrays
#define VECTOR3R_TO_NUMPY(vec,arr) arr[0]=vec[0]; arr[1]=vec[1]; arr[2]=vec[2]
-
-#define MATRIX3R_TO_NUMPY(mat,arr) arr[0]=mat(0,0);arr[1]=mat(0,1);arr[2]=mat(0,2);arr[3]=mat(1,0);arr[4]=mat(1,1);arr[5]=mat(1,2);arr[6]=mat(2,0);arr[7]=mat(2,1);arr[8]=mat(2,2);
\ No newline at end of file
+#define MATRIX3R_TO_NUMPY(mat,arr) arr[0]=mat(0,0);arr[1]=mat(0,1);arr[2]=mat(0,2);arr[3]=mat(1,0);arr[4]=mat(1,1);arr[5]=mat(1,2);arr[6]=mat(2,0);arr[7]=mat(2,1);arr[8]=mat(2,2);
=== modified file 'pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.hpp'
--- pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.hpp 2010-01-11 12:05:56 +0000
+++ pkg/common/Engine/GlobalEngine/PersistentTriangulationCollider.hpp 2010-02-18 22:47:05 +0000
@@ -7,7 +7,6 @@
*************************************************************************/
#pragma once
-
#include<yade/core/Collider.hpp>
#include<yade/core/InteractionContainer.hpp>
#include <list>
@@ -15,9 +14,6 @@
#include <vector>
#include <algorithm>
-#include<yade/pkg-dem/TesselationWrapper.hpp>
-
-
/*! \brief Collision detection engine based on regular triangulation.
This engine is using CGAL library (see http://www.cgal.org/)
@@ -26,8 +22,8 @@
Also needed : uncommenting lines in core/Interaction.cpp and core/Interaction.hpp (see NOTE:TriangulationCollider
*/
+class TesselationWrapper;
-using namespace std;
class PersistentTriangulationCollider : public Collider
{
@@ -59,22 +55,22 @@
unsigned int nbObjects;
/// Aabb extremity of the sphere number "id" projected onto the X axis
- vector<shared_ptr<AABBBound> > xBounds;
+ std::vector<shared_ptr<AABBBound> > xBounds;
/// Aabb extremity of the sphere number "id" projected onto the Y axis
- vector<shared_ptr<AABBBound> > yBounds;
+ std::vector<shared_ptr<AABBBound> > yBounds;
/// Aabb extremity of the sphere number "id" projected onto the Z axis
- vector<shared_ptr<AABBBound> > zBounds;
+ std::vector<shared_ptr<AABBBound> > zBounds;
// collection of Aabb that are in interaction
//protected : vector< set<unsigned int> > overlappingBB;
shared_ptr<InteractionContainer> interactions;
/// upper right corner of the Aabb of the objects => for spheres = center[i]-radius
- vector<Real> maxima;
+ std::vector<Real> maxima;
/// lower left corner of the Aabb of the objects => for spheres = center[i]+radius
- vector<Real> minima;
+ std::vector<Real> minima;
/// Used the first time broadInteractionTest is called, to initialize and sort the xBounds, yBounds,
/// and zBounds arrays and to initialize the overlappingBB collection
=== modified file 'pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp'
--- pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp 2010-02-14 21:23:14 +0000
+++ pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.cpp 2010-02-18 22:47:05 +0000
@@ -10,30 +10,13 @@
#include<yade/pkg-dem/FrictPhys.hpp>
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-dem/DemXDofGeom.hpp>
-//#include<yade/pkg-dem/MacroMicroElasticRelationships.hpp>
#include<yade/core/Interaction.hpp>
#include<yade/core/Scene.hpp>
CREATE_LOGGER(GlobalStiffnessTimeStepper);
YADE_PLUGIN((GlobalStiffnessTimeStepper));
-GlobalStiffnessTimeStepper::GlobalStiffnessTimeStepper() : TimeStepper()
-{
-//cerr << "GlobalStiffnessTimeStepper()" << endl;
- sdecGroupMask = 1;
- timestepSafetyCoefficient = 0.8;
- computedOnce = false;
- defaultDt = 1;
- previousDt = defaultDt;
-
-}
-
-
-GlobalStiffnessTimeStepper::~GlobalStiffnessTimeStepper()
-{
-
-}
-
-
+
+GlobalStiffnessTimeStepper::~GlobalStiffnessTimeStepper() {}
void GlobalStiffnessTimeStepper::findTimeStepFromBody(const shared_ptr<Body>& body, Scene * ncb)
{
@@ -41,28 +24,10 @@
Vector3r& stiffness= stiffnesses[body->getId()];
Vector3r& Rstiffness=Rstiffnesses[body->getId()];
-
- if(! ( /* sphere && */ sdec && stiffness==Vector3r::ZERO) )
+
+ if(!sdec || stiffness==Vector3r::ZERO)
return; // not possible to compute!
- //cerr << "return; // not possible to compute!" << endl;
-// Real Dab = sphere->radius;
-// Real rad3 = std::pow(Dab,2); // radius to the power of 2, from sphere
-//
-// Real Eab = sdec->young;
-// Real Vab = sdec->poisson;
-// Real Dinit = 2*Dab; // assuming thet sphere is in contact with itself
-// Real Sinit = Mathr::PI * std::pow( Dab , 2);
-//
-// Real alpha = sdecContactModel->alpha;
-// Real beta = sdecContactModel->beta;
-// Real gamma = sdecContactModel->gamma;
-//
-// Real Kn = abs((Eab*Sinit/Dinit)*( (1+alpha)/(beta*(1+Vab) + gamma*(1-alpha*Dab) ) ));
-// Real Ks = abs(Kn*(1-alpha*Vab)/(1+Vab));
-
-// Real dt = min (sqrt( sdec->mass / stiffness.x()) ,sqrt( sdec->mass / stiffness.y() ) );
-// dt = 0.5 * min (sqrt( sdec->mass / stiffness.z()), dt);
Real dtx, dty, dtz;
Real dt = max( max (stiffness.X(), stiffness.Y()), stiffness.Z() );
@@ -73,7 +38,7 @@
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;
@@ -81,30 +46,10 @@
if (Rstiffness.Z()!=0) {
dtz = sdec->inertia.Z()/Rstiffness.Z(); computedSomething = true;}
else dtz = Mathr::MAX_REAL;
-
-
- //if (Rstiffness.X()!=0)
- //Real dtx = (Rstiffness.X()!=0 ? sdec->inertia.X()/Rstiffness.X() : Mathr::MAX_REAL);
- //Real dty = (Rstiffness.Y()!=0 ? sdec->inertia.Y()/Rstiffness.Y() : Mathr::MAX_REAL);
- //Real dtz = (Rstiffness.Z()!=0 ? sdec->inertia.Z()/Rstiffness.Z() : Mathr::MAX_REAL);
+
Real Rdt = std::min( std::min (dtx, dty), dtz );//Rdt = smallest squared eigenperiod for rotational motions
-
- //cerr << "Rstiffness.x()=" << Rstiffness.x() << " " << Rstiffness.y() << " " << Rstiffness.z() << endl;
- //cerr << "sdec->inertia=" << sdec->inertia.x() << " " << sdec->inertia.x() << " " << sdec->inertia.x() << endl;
- //cerr << "timesteps : dt=" << dt << " / Rdt=" << Rdt << endl;
-
dt = 1.41044*timestepSafetyCoefficient*std::sqrt(std::min(dt,Rdt));//1.41044 = sqrt(2)
-
newDt = std::min(dt,newDt);
- //computedSomething = true;
-
-}
-
-
-void GlobalStiffnessTimeStepper::findTimeStepFromInteraction(const shared_ptr<Interaction>& interaction, shared_ptr<BodyContainer>& bodies)
-{
-
-
}
bool GlobalStiffnessTimeStepper::isActivated(Scene*)
@@ -112,11 +57,11 @@
return (active && ((!computedOnce) || (Omega::instance().getCurrentIteration() % timeStepUpdateInterval == 0) || (Omega::instance().getCurrentIteration() < (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.
+ // ANSW : Really strange!! I also noticed timestepper was not computing anything if I don't set active=true in python, why? since the default is true, it seems. I didn't try the compare script yet. (BC.)
active=active;
computeStiffnesses(ncb);
@@ -124,44 +69,21 @@
// shared_ptr<InteractionContainer>& interactions = ncb->interactions;
newDt = Mathr::MAX_REAL;
- //Real defaultDt = 0.0003;
-
- //computedSomething = false; // this flag is to avoid setting timestep to MAX_REAL :)
-/*
- InteractionContainer::iterator ii = interactions->begin();
- InteractionContainer::iterator iiEnd = interactions->end();
- for( ; ii!=iiEnd ; ++ii )
- findTimeStepFromInteraction(*ii , bodies);*/
-
- //if(!computedSomething)
- //{
- BodyContainer::iterator bi = bodies->begin();
- BodyContainer::iterator biEnd = bodies->end();
- for( ; bi!=biEnd ; ++bi )
- {
- //cerr << "for( ; bi!=biEnd ; ++bi )" << endl;
- shared_ptr<Body> b = *bi;
- if (b->isDynamic) {
- //cerr << "if (body->isDynamic) {" << endl;
- if((b->getGroupMask() & sdecGroupMask)||sdecGroupMask==0)
- //cerr << "if( b->getGroupMask() & sdecGroupMask)" << computedSomething << endl;
- findTimeStepFromBody(b, ncb); }
- }
- //}
+ BodyContainer::iterator bi = bodies->begin();
+ BodyContainer::iterator biEnd = bodies->end();
+ 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
Omega::instance().setTimeStep(previousDt);
- //previousDt = Omega::instance().getTimeStep();
- //Omega::instance().setTimeStep(newDt);
computedOnce = true;
- //cerr << "computedOnce=" << computedOnce << endl;
- //cerr << "computed timestep is:" << newDt;
}
else if (!computedOnce) Omega::instance().setTimeStep(defaultDt);
- //cerr << " new timestep is:" << Omega::instance().getTimeStep() << endl;
-
LOG_INFO("computed timestep " << newDt <<
(Omega::instance().getTimeStep()==newDt ? string(", appplied") :
string(", BUT timestep is ")+lexical_cast<string>(Omega::instance().getTimeStep()))<<".");
@@ -184,33 +106,26 @@
NormShearPhys* phys=YADE_CAST<NormShearPhys*>(contact->interactionPhysics.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;
- // FIXME? NormShearPhys knows nothing about whether the contact is "active" (force!=0) or not;
- // former code: if(force==0) continue; /* disregard this interaction, it is not active */.
- // It seems though that in such case either the interaction is accidentally at perfect equilibrium (unlikely)
- // or it should have been deleted already. Right?
- //ANSWER : some interactions can exist without fn, e.g. distant capillary force, wich does not contribute to the overall stiffness via kn. The test is needed.
Real fn = (static_cast<NormShearPhys *> (contact->interactionPhysics.get()))->normalForce.SquaredLength();
-
- if (fn!=0) {
- //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;
- //cerr << "diag_Rstifness=" << diag_Rstiffness << endl;
-
- 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);
- }
+ 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;
+ //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);
}
}
=== modified file 'pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp'
--- pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp 2009-12-04 23:07:34 +0000
+++ pkg/dem/Engine/GlobalEngine/GlobalStiffnessTimeStepper.hpp 2010-02-18 22:47:05 +0000
@@ -12,14 +12,11 @@
/*! \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)
+ See usage details in TriaxialTest documentation (TriaxialTest is also a good example of how to use this class)
*/
class Interaction;
class BodyContainer;
-//class MacroMicroElasticRelationships;
class Scene;
class GlobalStiffnessTimeStepper : public TimeStepper
@@ -29,31 +26,25 @@
vector<Vector3r> Rstiffnesses;
void computeStiffnesses(Scene*);
- Real newDt, previousDt;
+ Real newDt;
bool computedSomething,
- computedOnce;
- //shared_ptr<MacroMicroElasticRelationships> sdecContactModel;
-
+ computedOnce;
void findTimeStepFromBody(const shared_ptr<Body>& body, Scene * ncb);
- void findTimeStepFromInteraction(const shared_ptr<Interaction>& , shared_ptr<BodyContainer>&);
-
+
public :
- int sdecGroupMask; // FIXME - we should find a way to clean groupmask stuff
- //! defaultDt is used as default AND as max value of the timestep
- Real defaultDt;
+ int sdecGroupMask; // FIXME - to be removed -> not used here but set in preprocessors, removing breaks compilation.
//! used as a multiplier on the theoretical critical timestep (compensate some approximations in the computation)
- Real timestepSafetyCoefficient;
-
- GlobalStiffnessTimeStepper();
virtual ~GlobalStiffnessTimeStepper();
virtual void computeTimeStep(Scene*);
virtual bool isActivated(Scene*);
+ 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,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;)
DECLARE_LOGGER;
-
- REGISTER_ATTRIBUTES(TimeStepper,(sdecGroupMask)(defaultDt)(previousDt)(timestepSafetyCoefficient)(computedOnce));
- REGISTER_CLASS_NAME(GlobalStiffnessTimeStepper);
- REGISTER_BASE_CLASS_NAME(TimeStepper);
};
REGISTER_SERIALIZABLE(GlobalStiffnessTimeStepper);
=== modified file 'pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp'
--- pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp 2010-02-15 17:51:06 +0000
+++ pkg/dem/Engine/GlobalEngine/TesselationWrapper.cpp 2010-02-18 22:47:05 +0000
@@ -102,7 +102,7 @@
else {
v->info() = (const unsigned int) p->second;
//Vh->info().isFictious = false;//false is the default
- Tes.max_id = std::max(Tes.max_id, p->second);
+ Tes.max_id = std::max(Tes.max_id,(int) p->second);
Tes.vertexHandles[p->second]=v;
hint=v->cell();
++TW.n_spheres;
