yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #06072
[Branch ~yade-dev/yade/trunk] Rev 2528: Re-write of "NormalInelasticityLaw" with use of ScGeom6D for what concerns rotations.
------------------------------------------------------------
revno: 2528
committer: jduriez <jduriez@c1solimara-l>
branch nick: yade
timestamp: Wed 2010-11-03 16:50:02 +0100
message:
Re-write of "NormalInelasticityLaw" with use of ScGeom6D for what concerns rotations.
=> Rename of the corresponding Law2...
Add a reference of a article. In references.bib AND in yade-articles.bib, because I understood that "references" is for what is quoted from the documentation, and "yade-articles" is for what was made with Yade. And here it concerns both...
Correct of two typos in doc of StepDisplacer
modified:
doc/references.bib
doc/yade-articles.bib
pkg/common/StepDisplacer.hpp
pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp
pkg/dem/KinemCNLEngine.hpp
pkg/dem/KinemCNSEngine.cpp
pkg/dem/KinemSimpleShearBox.cpp
pkg/dem/NormalInelasticityLaw.cpp
pkg/dem/NormalInelasticityLaw.hpp
pkg/dem/NormalInelasticityPhys.hpp
pkg/dem/SimpleShear.cpp
py/system.py
scripts/normalInelasticity-test.py
scripts/simpleShear.py
--
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 'doc/references.bib'
--- doc/references.bib 2010-10-13 15:11:01 +0000
+++ doc/references.bib 2010-11-03 15:50:02 +0000
@@ -37,6 +37,14 @@
doi = "10.1061/(ASCE)0733-9399(2005)131:7(689)"
}
+@Article{Duriez2010,
+ author = "J. Duriez and F.Darve and F.-V.Donze",
+ title = "A discrete modeling-based constitutive relation for infilled rock joints",
+ year = "2010",
+ journal = "International Journal of Rock Mechanics & Mining Sciences",
+ note = "in press"
+}
+
@PhDThesis{ Chareyre2003,
author = "Bruno Chareyre",
title = "Mod{\'e}lisation du comportement d'ouvrages composites sol-g{\'e}osynth{\'e}tique par {\'e}l{\'e}ments discrets - Application aux tranch{\'e}es d'ancrage en t{\^e}te de talus.",
=== modified file 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2010-10-20 11:31:49 +0000
+++ doc/yade-articles.bib 2010-11-03 15:50:02 +0000
@@ -271,3 +271,12 @@
journal = "Granular Matter",
year = "2010 (submitted)"
}
+
+@Article{Duriez2010,
+ author = "J. Duriez and F.Darve and F.-V.Donze",
+ title = "A discrete modeling-based constitutive relation for infilled rock joints",
+ year = "2010",
+ journal = "International Journal of Rock Mechanics & Mining Sciences",
+ note = "in press"
+}
+
=== modified file 'pkg/common/StepDisplacer.hpp'
--- pkg/common/StepDisplacer.hpp 2010-10-13 16:23:08 +0000
+++ pkg/common/StepDisplacer.hpp 2010-11-03 15:50:02 +0000
@@ -7,7 +7,7 @@
public:
virtual void action();
YADE_CLASS_BASE_DOC_ATTRS(StepDisplacer,PartialEngine,"Apply generalized displacement (displacement or rotation) stepwise on subscribed bodies.",
- ((Se3r,deltaSe3,Se3r(Vector3r(NaN,NaN,NaN),Quaternionr::Identity()),,"|ydeprec| Set moc/rot directly instead; keps only for backwards compat. If the 0th component of the vector is not NaN, then it was updated by the user, warning is issued and mov and rot are updated automatically."))
+ ((Se3r,deltaSe3,Se3r(Vector3r(NaN,NaN,NaN),Quaternionr::Identity()),,"|ydeprec| Set mov/rot directly instead; kept only for backwards compat. If the 0th component of the vector is not NaN, then it was updated by the user, warning is issued and mov and rot are updated automatically."))
((Vector3r,mov,Vector3r::Zero(),,"Linear displacement step to be applied per iteration, by addition to :yref:`State.pos`."))
((Quaternionr,rot,Quaternionr::Identity(),,"Rotation step to be applied per iteration (via rotation composition with :yref:`State.ori`)."))
((bool,setVelocities,false,,"If true, velocity and angularVelocity are modified in such a way that over one iteration (dt), the body will have the prescribed jump. In this case, :yref:`position<State.pos>` and :yref:`orientation<State.ori>` itself is not updated for :yref:`dynamic<Body::dynamic>` bodies, since they would have the delta applied twice (here and in the :yref:`integrator<NewtonIntegrator>`). For non-dynamic bodies however, they *are* still updated."))
=== modified file 'pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp'
--- pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp 2010-11-03 15:50:02 +0000
@@ -10,7 +10,7 @@
#include<yade/pkg-common/Dispatching.hpp>
-/*! \brief The RelationShips for using Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity
+/*! \brief The RelationShips for using Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity
In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed.
WARNING : as in the others Relationships most of the attributes are computed only once : when the interaction is "new"
@@ -27,7 +27,7 @@
FUNCTOR2D(NormalInelasticMat,NormalInelasticMat);
YADE_CLASS_BASE_DOC_ATTRS(Ip2_2xNormalInelasticMat_NormalInelasticityPhys,
IPhysFunctor,
- "The RelationShips for using Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity\n\n In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
+ "The RelationShips for using Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity\n\n In these RelationShips all the attributes of the interactions (which are of NormalInelasticityPhys type) are computed. \n\n.. warning::\n\tas in the others :yref:`Ip2 functors<IPhysFunctor>`, most of the attributes are computed only once, when the interaction is new.",
((Real,betaR,0.12,,"Parameter for computing the torque-stifness : T-stifness = betaR * Rmoy^2"))
);
};
=== modified file 'pkg/dem/KinemCNLEngine.hpp'
--- pkg/dem/KinemCNLEngine.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/KinemCNLEngine.hpp 2010-11-03 15:50:02 +0000
@@ -32,7 +32,7 @@
((Real,gammalim,0.0,,"the value of tangential displacement (of upper plate) at wich the shearing is stopped [m]"))
((Real,gamma,0.0,,"current value of tangential displacement [m]"))
((std::vector<Real>,gamma_save,,,"vector with the values of gamma at which a save of the simulation is performed [m]")),
-// ((Real,coeff_dech,1.0,,"in the case of the use of 'Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity' for ex, where kn(unload)#kn(load). The engine cares to find the value at the first run BROKEN actually")),
+// ((Real,coeff_dech,1.0,,"in the case of the use of 'Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity' for ex, where kn(unload)#kn(load). The engine cares to find the value at the first run BROKEN actually")),
temoin=0;
it_stop=0;
);
=== modified file 'pkg/dem/KinemCNSEngine.cpp'
--- pkg/dem/KinemCNSEngine.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/KinemCNSEngine.cpp 2010-11-03 15:50:02 +0000
@@ -111,11 +111,11 @@
// // vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
// // for ( ;itFirst!=itLast; ++itFirst )
// // {
-// // if ( ( *itFirst )->getClassName() == "Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity" )
+// // if ( ( *itFirst )->getClassName() == "Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity" )
// // {
-// // myLdc = YADE_PTR_CAST<Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity> ( *itFirst );
+// // myLdc = YADE_PTR_CAST<Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity> ( *itFirst );
// // coeff_dech = myLdc ->coeff_dech;
-// // if(LOG) cout << "Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity engine found, with coeff_dech = " << coeff_dech << endl;
+// // if(LOG) cout << "Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity engine found, with coeff_dech = " << coeff_dech << endl;
// // }
// // }
// }
@@ -132,7 +132,7 @@
// Real Hcurrent = topbox->state->pos.y();
// Real Fdesired = F_0 + KnC * 1.0e9 * Scontact * (Hcurrent-Y0); // The value of the force desired
//
-// // Prise en compte de la difference de rigidite entre charge et decharge dans le cadre de Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity : => INUTILE maintenant ?
+// // Prise en compte de la difference de rigidite entre charge et decharge dans le cadre de Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : => INUTILE maintenant ?
// // if( F_sup.y() > Fdesired ) // cas ou l'on va monter la plaq <=> (normalemt) a une decharge
// // stiffness *= coeff_dech;
//
=== modified file 'pkg/dem/KinemSimpleShearBox.cpp'
--- pkg/dem/KinemSimpleShearBox.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/KinemSimpleShearBox.cpp 2010-11-03 15:50:02 +0000
@@ -169,11 +169,11 @@
// vector<shared_ptr<Engine> >::iterator itLast = scene->engines.end();
// for ( ;itFirst!=itLast; ++itFirst )
// {
-// if ( ( *itFirst )->getClassName() == "Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity" )
+// if ( ( *itFirst )->getClassName() == "Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity" )
// {
-// myLdc = YADE_PTR_CAST<Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity> ( *itFirst );
+// myLdc = YADE_PTR_CAST<Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity> ( *itFirst );
// coeff_dech = myLdc ->coeff_dech;
-// if(LOG) cout << "Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity engine found, with coeff_dech = " << coeff_dech << endl;
+// if(LOG) cout << "Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity engine found, with coeff_dech = " << coeff_dech << endl;
// }
// }
// }
@@ -194,7 +194,7 @@
// cout << "Alors que f0 = = "<< f0 << endl;
// cout << "Car terme correctif = " << KnC * 1.0e9 * Scontact * (hCurrent-y0)<< endl;
-// Prise en compte de la difference de rigidite entre charge et decharge dans le cadre de Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity : => INUTILE maintenant ?
+// Prise en compte de la difference de rigidite entre charge et decharge dans le cadre de Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : => INUTILE maintenant ?
// if( fSup.y() > fDesired ) // cas ou l'on va monter la plaq <=> (normalemt) a une decharge
// stiffness *= coeff_dech;
=== modified file 'pkg/dem/NormalInelasticityLaw.cpp'
--- pkg/dem/NormalInelasticityLaw.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NormalInelasticityLaw.cpp 2010-11-03 15:50:02 +0000
@@ -12,11 +12,11 @@
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-YADE_PLUGIN((Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity));
-
-
-
-void Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
+YADE_PLUGIN((Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity));
+
+
+
+void Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
{
const Real& dt = scene->dt;
@@ -28,7 +28,7 @@
State* de1 = Body::byId(id1,scene)->state.get();
State* de2 = Body::byId(id2,scene)->state.get();
NormalInelasticMat* Mat1 = static_cast<NormalInelasticMat*>(Body::byId(id1,scene)->material.get());
- ScGeom* currentContactGeometry = YADE_CAST<ScGeom*>(iG.get());
+ ScGeom6D* currentContactGeometry = YADE_CAST<ScGeom6D*>(iG.get());
NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (iP.get());
Vector3r& shearForce = currentContactPhysics->shearForce;
@@ -42,9 +42,7 @@
}
- // *** Computation of normal Force : depends of the history *** //
- if( currentContactGeometry->penetrationDepth < currentContactPhysics->unMax)
- currentContactPhysics->kn *= Mat1->coeff_dech;
+// ******** Computation of normal Force : depends of the history ******* //
Real Fn; // la valeur de Fn qui va etre calculee selon différentes manieres puis affectee
// cout << " Dans Law2 valeur de kn : " << currentContactPhysics->kn << endl;
@@ -58,12 +56,13 @@
currentContactPhysics->unMax = std::abs(un);
// cout << "je suis dans le calcul normal " << endl;
}
- else// a priori then we need a greater stifness. False in the case when we are already on the limit state, but a correction below will then do the jobthe law2 will do the job in the case with a correction
+ else// a priori then we need a greater stifness. False in the case when we are already on the limit state, but a correction below will then do the job
{
currentContactPhysics->kn = currentContactPhysics->knLower* Mat1->coeff_dech;
- Fn = currentContactPhysics->previousFn + currentContactPhysics->kn * (un-currentContactPhysics->previousun); // Calcul incrémental de la nvlle force
+ // Incremental computation of the new normal force :
+ Fn = currentContactPhysics->previousFn + currentContactPhysics->kn * (un-currentContactPhysics->previousun);
// cout << "je suis dans l'autre calcul" << endl;
- if(std::abs(Fn) > std::abs(currentContactPhysics->knLower * un)) // check the limit state is not violated
+ if(std::abs(Fn) > std::abs(currentContactPhysics->knLower * un)) // check if the limit state is not violated
{
Fn = currentContactPhysics->knLower*un;
// cout << "j'etais dans l'autre calcul mais j'ai corrige pour ne pas depasser la limite" << endl;
@@ -83,7 +82,9 @@
// *** End of computation of normal force *** //
-
+
+
+// ******** Tangential force ******* //
if ( un < 0 )
{ // BREAK due to tension
scene->interactions->requestErase(contact->getId1(),contact->getId2());
@@ -119,15 +120,14 @@
// q.fromAngleAxis(angle,axis);
// currentContactPhysics->shearForce = q*currentContactPhysics->shearForce;
- /// ///
-
+
Vector3r x = currentContactGeometry->contactPoint;
Vector3r c1x = (x - de1->se3.position);
Vector3r c2x = (x - de2->se3.position);
/// The following definition of c1x and c2x is to avoid "granular ratcheting"
/// (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
/// "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
- /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors)
+ /// ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors, or discussions on yade mailing lists)
Vector3r _c1x_ = currentContactGeometry->radius1*currentContactGeometry->normal;
Vector3r _c2x_ = -currentContactGeometry->radius2*currentContactGeometry->normal;
Vector3r relativeVelocity = (de2->vel+de2->angVel.cross(_c2x_)) - (de1->vel+de1->angVel.cross(_c1x_));
@@ -143,12 +143,11 @@
if ( Fs > maxFs )
{
-
maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
maxFs = maxFs / Fs;
if (maxFs>1)
- cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
+ cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!" << endl;
shearForce *= maxFs;
if (Fn<0) currentContactPhysics->normalForce = Vector3r::Zero();
}
@@ -159,37 +158,37 @@
scene->forces.addForce (id2, f);
scene->forces.addTorque(id1,-c1x.cross(f));
scene->forces.addTorque(id2, c2x.cross(f));
+// *** End of computation of tangential force *** //
// ******** Moment law ******* //
if(momentRotationLaw)
{
- {// updates only orientation of contact (local coordinate system)
- Vector3r axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal); axis.normalize();
- Real angle = unitVectorsAngle(currentContactPhysics->prevNormal,currentContactGeometry->normal);
- Quaternionr align(AngleAxisr(angle,axis));
- currentContactPhysics->currentContactOrientation = align * currentContactPhysics->currentContactOrientation;
- }
-
- Quaternionr delta( de1->se3.orientation * currentContactPhysics->initialOrientation1.conjugate() *
- currentContactPhysics->initialOrientation2 * de2->se3.orientation.conjugate());
-
- AngleAxisr aa(delta); // aa.axis() of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
- if(angle > Mathr::PI) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
+// {// updates only orientation of contact (local coordinate system)
+// Vector3r axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal); axis.normalize();
+// Real angle = unitVectorsAngle(currentContactPhysics->prevNormal,currentContactGeometry->normal);
+// Quaternionr align(AngleAxisr(angle,axis));
+// currentContactPhysics->currentContactOrientation = align * currentContactPhysics->currentContactOrientation;
+// }
+//
+// Quaternionr delta( de1->se3.orientation * currentContactPhysics->initialOrientation1.conjugate() *
+// currentContactPhysics->initialOrientation2 * de2->se3.orientation.conjugate());
+//
+// AngleAxisr aa(delta); // aa.axis() of rotation - this is the Moment direction UNIT vector; angle represents the power of resistant ELASTIC moment
+// if(angle > Mathr::PI) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
//Real elasticMoment = currentContactPhysics->kr * std::abs(angle); // positive value (*)
- Real angle_twist(aa.angle() * aa.axis().dot(currentContactGeometry->normal) );
- Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
- Vector3r moment_twist(axis_twist * currentContactPhysics->kr);
-
- Vector3r axis_bending(aa.angle()*aa.axis() - axis_twist);
- Vector3r moment_bending(axis_bending * currentContactPhysics->kr);
+// Real angle_twist(aa.angle() * aa.axis().dot(currentContactGeometry->normal) );
+// Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
+ Vector3r moment_twist( (currentContactGeometry->getTwist()*currentContactPhysics->kr)*currentContactGeometry->normal );
+// Vector3r axis_bending(aa.angle()*aa.axis() - axis_twist);
+ Vector3r moment_bending( currentContactGeometry->getBending() * currentContactPhysics->kr );
Vector3r moment = moment_twist + moment_bending;
-// Limitation par seuil plastique
+// Limitation by plastic threshold
if (!momentAlwaysElastic)
{
Real normeMomentMax = currentContactPhysics->forMaxMoment * std::fabs(Fn);
=== modified file 'pkg/dem/NormalInelasticityLaw.hpp'
--- pkg/dem/NormalInelasticityLaw.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NormalInelasticityLaw.hpp 2010-11-03 15:50:02 +0000
@@ -16,21 +16,21 @@
-class Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
+class Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
{
public :
virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
FUNCTOR2D(ScGeom,NormalInelasticityPhys);
- YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity,
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity,
LawFunctor,
- "Contact law used to simulate granulate filler in rock joints [Duriez2009a]_. It includes possibility of cohesion, moment transfer and inelastic compression behaviour (to reproduce the normal inelasticity observed for rock joints, for the latter).\n\n The moment transfer relation corresponds to the adaptation of the work of Plassiard & Belheine (see in [DeghmReport2006]_ for example), which was realized by J. Kozicki, and which is also present in other ''Law2'' (see http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg02444.html for a ''state of art'').\n\n As others :yref:`LawFunctor`, it uses pre-computed data of the interactions (rigidities, friction angles -with their tan()-, orientations of the interactions); this work is done here in :yref:`Ip2_2xNormalInelasticMat_NormalInelasticityPhys`.\n\n To use this you should also use :yref:`NormalInelasticMat` as material type of the bodies.\n\n The effects of this law are illustrated in scripts/normalInelasticityTest.py",
+ "Contact law used to simulate granulate filler in rock joints [Duriez2009a]_, [Duriez2010]_. It includes possibility of cohesion, moment transfer and inelastic compression behaviour (to reproduce the normal inelasticity observed for rock joints, for the latter).\n\n The moment transfer relation corresponds to the adaptation of the work of Plassiard & Belheine (see in [DeghmReport2006]_ for example), which was realized by J. Kozicki, and is now coded in ScGeom6D.\n\n As others :yref:`LawFunctor`, it uses pre-computed data of the interactions (rigidities, friction angles -with their tan()-, orientations of the interactions); this work is done here in :yref:`Ip2_2xNormalInelasticMat_NormalInelasticityPhys`.\n\n To use this you should also use :yref:`NormalInelasticMat` as material type of the bodies.\n\n The effects of this law are illustrated in scripts/normalInelasticityTest.py",
((bool,momentRotationLaw,true,,"boolean, true=> computation of a torque (against relative rotation) exchanged between particles"))
((bool,momentAlwaysElastic,false,,"boolean, true=> the torque (computed only if momentRotationLaw !!) is not limited by a plastic threshold"))
);
};
-REGISTER_SERIALIZABLE(Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity);
+REGISTER_SERIALIZABLE(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity);
=== modified file 'pkg/dem/NormalInelasticityPhys.hpp'
--- pkg/dem/NormalInelasticityPhys.hpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NormalInelasticityPhys.hpp 2010-11-03 15:50:02 +0000
@@ -17,7 +17,7 @@
virtual ~NormalInelasticityPhys();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
- "Physics (of interaction) for using :yref:`Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity` : with inelastic unloadings",
+ "Physics (of interaction) for using :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity` : with inelastic unloadings",
((Real,unMax,0.0,,"the maximum value of penetration depth of the history of this interaction"))
((Real,previousun,0.0,,"the value of this un at the last time step"))
((Real,previousFn,0.0,,"the value of the normal force at the last time step"))
=== modified file 'pkg/dem/SimpleShear.cpp'
--- pkg/dem/SimpleShear.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/SimpleShear.cpp 2010-11-03 15:50:02 +0000
@@ -220,7 +220,7 @@
scene->engines.push_back(collider);
scene->engines.push_back(interactionGeometryDispatcher);
scene->engines.push_back(interactionPhysicsDispatcher);
-// scene->engines.push_back(shared_ptr<Engine>(new Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity));
+// scene->engines.push_back(shared_ptr<Engine>(new Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity));
if(gravApplied)
scene->engines.push_back(gravityCondition);
scene->engines.push_back(shared_ptr<Engine> (new NewtonIntegrator));
=== modified file 'py/system.py'
--- py/system.py 2010-10-30 17:08:51 +0000
+++ py/system.py 2010-11-03 15:50:02 +0000
@@ -145,6 +145,7 @@
'Law2_ScGeom_FrictPhys_Basic':'Law2_ScGeom_FrictPhys_CundallStrack', # Wed Oct 13 17:40:42 2010, bchareyre@dt-rv020
'Law2_Dem3DofGeom_FrictPhys_Basic':'Law2_Dem3DofGeom_FrictPhys_CundallStrack', # Wed Oct 13 17:45:32 2010, bchareyre@dt-rv020
'Law2_ScGeom_CohFrictPhys_ElasticPlastic':'Law2_ScGeom_CohFrictPhys_CohesionMoment', # Wed Oct 13 17:47:09 2010, bchareyre@dt-rv020
+ 'Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity':'Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity', # Thu Oct 28 18:09:16 2010, jduriez@c1solimara-l
'SpiralEngine':'HelixEngine', # Sat Oct 30 05:52:06 2010, vaclav@flux
'InterpolatingSpiralEngine':'InterpolatingHelixEngine', # Sat Oct 30 05:52:21 2010, vaclav@flux
'SpiralInteractionLocator2d':'HelixInteractionLocator2d', # Sat Oct 30 05:53:04 2010, vaclav@flux
=== modified file 'scripts/normalInelasticity-test.py'
--- scripts/normalInelasticity-test.py 2010-09-27 17:47:59 +0000
+++ scripts/normalInelasticity-test.py 2010-11-03 15:50:02 +0000
@@ -5,7 +5,10 @@
#- then, test in tangential direction
#- finally for what concerns moment transfer, with relative orientation changes (use frictionAngle=0.0 in this case to study more easily what's going on)
-#Different graphs illustrate the effects of the different loadings. The run is paused at each plot window (so that there is time to observe them). Type "Return" in the Yade terminal to go ahead.
+#Different graphs illustrate the effects of the different loadings. The run is paused at each plot window (so that there is time to observe them). Push on "Return", while being in the Yade terminal, to go ahead.
+
+#All comments concern r2528
+
from yade import plot
@@ -27,7 +30,7 @@
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_2xNormalInelasticMat_NormalInelasticityPhys(betaR=0.24)],
- [Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity()]
+ [Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity()]
),
PyRunner(iterPeriod=1,command='letMove()')
]
@@ -55,7 +58,7 @@
# ------ Test of the law in the normal direction, using python commands to let move ------ #
-
+print 'Beginning of normal loading'
O.dt=1e-5
yade.qt.View()
@@ -70,7 +73,8 @@
plot.plots={'step':('unVrai',),'unPerso':('normFn',),'unVrai':('normFnBis',)}
plot.plot()
raw_input()
-
+print 'End of normal loading'
+print ''
#NB : these different unVrai and unPerso illustrate the definition of penetrationDepth really used in the code (computed in Ig2_Sphere_Sphere_ScGeom) which is slightly different from R1 + R2 - Distance (see for example this "shift2"). According to the really used penetrationDepth, Fn evolves as it should
#O.saveTmp('EndComp')
@@ -79,10 +83,11 @@
# ------ Test of the law in the tangential direction, using StepDisplacer ------ #
+print 'Beginning of tangential loading'
+
dpos=Vector3.Zero
Vector3.__init__(dpos,1*O.dt,0,0)
-
-O.engines=O.engines[:4]+[StepDisplacer(subscribedBodies=[1],deltaSe3=(dpos,Quaternion.Identity),setVelocities=True)]+O.engines[5:]
+O.engines=O.engines[:3]+[StepDisplacer(ids=[1],mov=dpos,setVelocities=True)]+O.engines[4:]
O.run(1000)
plot.plots={'step':('gamma',),'gamma':('fx',)}
plot.plot()
@@ -92,19 +97,20 @@
raw_input()
#pylab.show() #to pause on the plot window. Effective only first time
-#-- Comments --#
+print 'End of tangential loading'
+print ''
+
+#-- Comments (r2528) --#
# - evolution of Fx with gamma normal (flat at the beginning because of the order of engines)
# - un decreases indeed during this shear, but maybe a zoom on the curves is needed to see it.
# - We can observe that the force state of the sample decreases a line with a slope equal to tan(~34.5°)=tan(~0.602 rad). Why not strict equality ? Because of the measure of the slope or because something else ? To see...
# - during this phase O.forces.t(0).norm() / O.forces.f(0)[0] seems to increase between 0.502 and 0.507 (according to r=[T[i]/F[i] for i in range(50,T.__len__()-20) ])
-# Note to explain this that Fx = O.forces.f(0)[0] is more and more different from Ft, from which we can expect Ft =
-# Torque /Radius
-
-
+# Note to explain this that Fx = O.forces.f(0)[0] is more and more different from Ft, from which we can expect Ft = Torque /Radius
## ------ Test of the law for the moment, using blockedDOF_s ------ #
#O.loadTmp('EndComp')
+print 'Beginning of rotationnal loading'
##To use blockedDOF_s, the body has to be dynamic....
upperSphere.dynamic=True
upperSphere.state.blockedDOFs='x','rx','y','ry','z','rz'
@@ -119,7 +125,8 @@
print upperSphere.state.ori
print i.geom.penetrationDepth
-
O.run(8000,True)
plot.plots={'step':('torque',)}
plot.plot()
+
+#-- Comments : TO DO
=== modified file 'scripts/simpleShear.py'
--- scripts/simpleShear.py 2010-09-27 17:47:59 +0000
+++ scripts/simpleShear.py 2010-11-03 15:50:02 +0000
@@ -65,7 +65,7 @@
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_2xNormalInelasticMat_NormalInelasticityPhys()],
- [Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity()]
+ [Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity()]
),
NewtonIntegrator(damping=.2)
,PyRunner(iterPeriod=50,command='defData()')
