yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04775
[Branch ~yade-dev/yade/trunk] Rev 2277: - Cleaning changes in NormalInelasticityFiles...
------------------------------------------------------------
revno: 2277
committer: jduriez <jduriez@c1solimara-l>
branch nick: trunk
timestamp: Fri 2010-06-04 17:26:30 +0200
message:
- Cleaning changes in NormalInelasticityFiles...
=>Add of a new Material storing the info for differences between loading and unloading (related changes)
- Suppress of initialKn/s in FrictPhys, (except in in attic...)
added:
pkg/dem/DataClass/Material/NormalInelasticMat.cpp
pkg/dem/DataClass/Material/NormalInelasticMat.hpp
renamed:
pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.cpp => pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp
pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp => pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp
modified:
pkg/dem/DataClass/InteractionPhysics/FrictPhys.hpp
pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.cpp
pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp
pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp
pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp
pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp
pkg/dem/PreProcessor/SimpleShear.cpp
py/system.py
scripts/NormalInelasticityTest.py
pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp
pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.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/DataClass/InteractionPhysics/FrictPhys.hpp'
--- pkg/dem/DataClass/InteractionPhysics/FrictPhys.hpp 2010-04-26 09:02:51 +0000
+++ pkg/dem/DataClass/InteractionPhysics/FrictPhys.hpp 2010-06-04 15:26:30 +0000
@@ -13,12 +13,10 @@
{
public :
// kn,ks,normal inherited from NormShearPhys
- Real initialKn // initial normal elastic constant.
- ,initialKs // initial shear elastic constant.
+ Real frictionAngle // angle of friction, according to Coulumb criterion
//Those two attributes were removed. If you need them, copy those two lines in your class inheriting from FrictPhys
// ,equilibriumDistance // equilibrium distance
// ,initialEquilibriumDistance // initial equilibrium distance
- ,frictionAngle // angle of friction, according to Coulumb criterion
;
virtual ~FrictPhys();
YADE_CLASS_BASE_DOC_ATTRS_CTOR(FrictPhys,NormShearPhys,"Interaction with friction",
=== modified file 'pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.cpp'
--- pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.cpp 2010-02-26 23:29:53 +0000
+++ pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.cpp 2010-06-04 15:26:30 +0000
@@ -9,15 +9,6 @@
#include "NormalInelasticityPhys.hpp"
-
-void NormalInelasticityPhys::SetBreakingState()
-{
-
- //if (fragile) {
-
-
-}
-
NormalInelasticityPhys::~NormalInelasticityPhys()
{
}
=== modified file 'pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp'
--- pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp 2010-06-03 10:12:50 +0000
+++ pkg/dem/DataClass/InteractionPhysics/NormalInelasticityPhys.hpp 2010-06-04 15:26:30 +0000
@@ -8,23 +8,16 @@
#pragma once
-#include<yade/pkg-dem/FrictPhys.hpp> // ou
-// #include <yade/pkg/dem/DataClass/InteractionPhysics/FrictPhys.hpp>
-
-/*! \brief Interaction for using Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity
-
-This interaction is similar to CohFrictPhys. Among the differences are the unMax, previousun and previousFn (allowing to describe the inelastic unloadings in compression), no more shear and tension Adhesion, no more "fragile", "cohesionBroken" and "cohesionDisablesFriction"
- */
+#include<yade/pkg-dem/FrictPhys.hpp>
+
class NormalInelasticityPhys : public FrictPhys
{
public :
-
-// NormalInelasticityPhys();
virtual ~NormalInelasticityPhys();
- void SetBreakingState ();
+
YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
- "Physics (of interaction) for using :yref:`Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity`",
+ "Physics (of interaction) for using :yref:`Law2_ScGeom_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"))
@@ -37,7 +30,8 @@
((Vector3r,initialPosition1,Vector3r::Zero(),""))
((Vector3r,initialPosition2,Vector3r::Zero(),""))
((Real,forMaxMoment,1.0,"parameter stored for each interaction, and allowing to compute the maximum value of the exchanged torque : TorqueMax= forMaxMoment * NormalForce"))
- ((Real,kr,0.0,"the rolling stiffness of the rigidity")),
+ ((Real,kr,0.0,"the rolling stiffness of the interaction"))
+ ((Real,knLower,0.0,"the stifness corresponding to a virgin load for example")),
createIndex();
);
REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhys);
=== added file 'pkg/dem/DataClass/Material/NormalInelasticMat.cpp'
--- pkg/dem/DataClass/Material/NormalInelasticMat.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/DataClass/Material/NormalInelasticMat.cpp 2010-06-04 15:26:30 +0000
@@ -0,0 +1,12 @@
+
+#include "NormalInelasticMat.hpp"
+
+
+NormalInelasticMat::~NormalInelasticMat()
+{
+}
+
+YADE_PLUGIN((NormalInelasticMat));
+
+
+
=== added file 'pkg/dem/DataClass/Material/NormalInelasticMat.hpp'
--- pkg/dem/DataClass/Material/NormalInelasticMat.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/DataClass/Material/NormalInelasticMat.hpp 2010-06-04 15:26:30 +0000
@@ -0,0 +1,32 @@
+/*************************************************************************
+* Copyright (C) 2010 by Jerome Duriez <jerome.duriez@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/pkg-common/ElastMat.hpp>
+
+
+class NormalInelasticMat : public FrictMat
+{
+ public :
+ virtual ~NormalInelasticMat ();
+
+/// Serialization
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticMat,FrictMat,"Material class for particles whose contact obey to a normal inelasticity (governed by this *coeff_dech*).",
+ ((Real,coeff_dech,1.0,"=kn(unload) / kn(load)"))
+ ,
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(NormalInelasticMat,FrictMat);
+};
+
+REGISTER_SERIALIZABLE(NormalInelasticMat);
+
+
=== renamed file 'pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.cpp' => 'pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 2010-06-04 15:26:30 +0000
@@ -6,10 +6,10 @@
* GNU General Public License v2 or later. See file LICENSE for details. *
*************************************************************************/
-#include"Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp"
+#include"Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp"
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-dem/NormalInelasticityPhys.hpp>
-#include<yade/pkg-dem/CohFrictMat.hpp>
+#include<yade/pkg-dem/NormalInelasticMat.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
@@ -23,21 +23,14 @@
//
//
-void Ip2_2xCohFrictMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // CohFrictMat
- , const shared_ptr<Material>& b2 // CohFrictMat
+void Ip2_2xNormalInelasticMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // NormalInelasticMat
+ , const shared_ptr<Material>& b2 // NormalInelasticMat
, const shared_ptr<Interaction>& interaction)
{
- CohFrictMat* sdec1 = static_cast<CohFrictMat*>(b1.get());
- CohFrictMat* sdec2 = static_cast<CohFrictMat*>(b2.get());
+ NormalInelasticMat* sdec1 = static_cast<NormalInelasticMat*>(b1.get());
+ NormalInelasticMat* sdec2 = static_cast<NormalInelasticMat*>(b2.get());
ScGeom* interactionGeometry = YADE_CAST<ScGeom*>(interaction->interactionGeometry.get());
- //Create cohesive interractions only once
- if (setCohesionNow && cohesionDefinitionIteration==-1) {
- cohesionDefinitionIteration=Omega::instance().getCurrentIteration();}
- if (setCohesionNow && cohesionDefinitionIteration!=-1 && cohesionDefinitionIteration!=Omega::instance().getCurrentIteration()) {
- cohesionDefinitionIteration = -1;
- setCohesionNow = 0;}
-
if(interactionGeometry) // so it is ScGeom - NON PERMANENT LINK
{
@@ -57,50 +50,26 @@
Real fb = sdec2->frictionAngle;
Real Kn = 2.0*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
+
Real Ks = 2.0*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
// Jean-Patrick Plassiard, Noura Belhaine, Frederic
- // Victor Donze, "Calibration procedure for spherical
+// Victor Donze, "Calibration procedure for spherical
// discrete elements using a local moemnt law".
Real Kr = betaR*std::pow((Da+Db)/2.0,2)*Ks;
- contactPhysics->initialKn = Kn;
- contactPhysics->initialKs = Ks;
-//cerr << "Ks: " << contactPhysics->initialKs << endl;
- contactPhysics->frictionAngle = std::min(fa,fb); // FIXME - this is actually a waste of memory space, just like initialKs and initialKn
+ contactPhysics->frictionAngle = std::min(fa,fb); // FIXME - this is actually a waste of memory space
contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
contactPhysics->forMaxMoment = 1.0*(Da+Db)/2.0; // 1.0 corresponding to ethaR which I don't know exactly where to define as a parameter...
- if ((setCohesionOnNewContacts || setCohesionNow) && sdec1->isCohesive && sdec2->isCohesive)
- {
+ // Lot of suppress here around (>) r2276. Normally not bad but ? See Ip2_2xCohFrictMat_CohFrictPhys.cpp to re-find the initial §...
- // FIXME - not sure: do I need to repeat it here [1] ?
- contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori;
- contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
- contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
- contactPhysics->initialPosition2 = Body::byId(interaction->getId2())->state->pos;
- contactPhysics->kr = Kr;
- contactPhysics->initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),interactionGeometry->normal);
- contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.conjugate() * contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.conjugate() * contactPhysics->initialContactOrientation;
-
-
-//if((interaction->getId1()==7 && interaction->getId2()==12)||(interaction->getId1()==12 && interaction->getId2()==7)){
-//Vector3r axis0;
-//Real angle0;
-//contactPhysics->initialRelativeOrientation.ToAxisAngle(axis0,angle0);
-//std::cout << "-----------------------\n"
-// << " ax0: " << axis0[0] << " " << axis0[1] << " " << axis0[2] << " an0: " << angle0 << "\n";
-//}
-
- //contactPhysics->elasticRollingLimit = elasticRollingLimit;
- }
-
contactPhysics->prevNormal = interactionGeometry->normal;
- contactPhysics->kn = contactPhysics->initialKn;
- contactPhysics->ks = contactPhysics->initialKs;
- // FIXME - or here [1] ?
+
+ contactPhysics->knLower = Kn;
+ contactPhysics->kn = Kn;
+
+ contactPhysics->ks = Ks;
contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori;
contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
@@ -110,56 +79,12 @@
contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.conjugate() * contactPhysics->initialContactOrientation;
contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.conjugate() * contactPhysics->initialContactOrientation;
- //contactPhysics->elasticRollingLimit = elasticRollingLimit;
- //
}
- else // !isNew
- {
- // FIXME - are those lines necessary ???? what they are doing in fact ???
- // ANSWER - they are used when you setCohesionNow (contact isNew not)
- NormalInelasticityPhys* contactPhysics = YADE_CAST<NormalInelasticityPhys*>(interaction->interactionPhysics.get());
-
- contactPhysics->kn = contactPhysics->initialKn;
- contactPhysics->ks = contactPhysics->initialKs;
- if (setCohesionNow && sdec1->isCohesive && sdec2->isCohesive)
- {
- //setCohesionNow = false;
-
- contactPhysics->initialOrientation1 = Body::byId(interaction->getId1())->state->ori;
- contactPhysics->initialOrientation2 = Body::byId(interaction->getId2())->state->ori;
- contactPhysics->initialPosition1 = Body::byId(interaction->getId1())->state->pos;
- contactPhysics->initialPosition2 = Body::byId(interaction->getId2())->state->pos;
- Real Da = interactionGeometry->radius1; // FIXME - multiply by factor of sphere interaction distance (so sphere interacts at bigger range that its geometrical size)
- Real Db = interactionGeometry->radius2; // FIXME - as above
- Real Kr = betaR*std::pow((Da+Db)/2.0,2)*contactPhysics->ks;
- contactPhysics->kr = Kr;
- contactPhysics->initialContactOrientation.setFromTwoVectors(Vector3r(1.0,0.0,0.0),interactionGeometry->normal);
- contactPhysics->currentContactOrientation = contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact1 = contactPhysics->initialOrientation1.conjugate() * contactPhysics->initialContactOrientation;
- contactPhysics->orientationToContact2 = contactPhysics->initialOrientation2.conjugate() * contactPhysics->initialContactOrientation;
-//Vector3r axis0;
-//Real angle0;
-//contactPhysics->initialRelativeOrientation.ToAxisAngle(axis0,angle0);
-//std::cout << "id1: " << interaction->getId1() << " id2: " << interaction->getId2() << " | "
-// << " ax0: " << axis0[0] << " " << axis0[1] << " " << axis0[2] << ", an0: " << angle0 << "\n";
- //contactPhysics->elasticRollingLimit = elasticRollingLimit;
- }
- }
}
-#if 0
- else // this is PERMANENT LINK because previous dynamic_cast failed, dispatcher should do this job
- {
- SDECLinkGeometry* sdecLinkGeometry = dynamic_cast<SDECLinkGeometry*>(interaction->interactionGeometry.get());
- if (sdecLinkGeometry)
- {
- SDECLinkPhysics* linkPhysics = static_cast<SDECLinkPhysics*>(interaction->interactionPhysics.get());
- // linkPhysics->frictionAngle = ?? //FIXME - uninitialized
- linkPhysics->kn = linkPhysics->initialKn;
- linkPhysics->ks = linkPhysics->initialKs;
- linkPhysics->equilibriumDistance = linkPhysics->initialEquilibriumDistance;
- }
- }
-#endif
+
};
-YADE_PLUGIN((Ip2_2xCohFrictMat_NormalInelasticityPhys));
+YADE_PLUGIN((Ip2_2xNormalInelasticMat_NormalInelasticityPhys));
+
+
+
=== renamed file 'pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp' => 'pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp'
--- pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp 2010-04-19 13:23:53 +0000
+++ pkg/dem/Engine/Functor/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp 2010-06-04 15:26:30 +0000
@@ -16,28 +16,22 @@
WARNING : as in the others Relationships most of the attributes are computed only once : when the interaction is "new"
*/
-class Ip2_2xCohFrictMat_NormalInelasticityPhys : public InteractionPhysicsFunctor
+class Ip2_2xNormalInelasticMat_NormalInelasticityPhys : public InteractionPhysicsFunctor
{
public :
-// Ip2_2xCohFrictMat_NormalInelasticityPhys();
virtual void go( const shared_ptr<Material>& b1,
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
-
- int cohesionDefinitionIteration; //useful is you want to use setCohesionNow
-
- FUNCTOR2D(CohFrictMat,CohFrictMat);
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Ip2_2xCohFrictMat_NormalInelasticityPhys,
+
+ FUNCTOR2D(NormalInelasticMat,NormalInelasticMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_2xNormalInelasticMat_NormalInelasticityPhys,
InteractionPhysicsFunctor,
"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<InteractionPhysicsFunctor>`, 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"))
- ((bool,setCohesionNow,false,""))
- ((bool,setCohesionOnNewContacts,false,"")),
- cohesionDefinitionIteration=-1;
);
};
-REGISTER_SERIALIZABLE(Ip2_2xCohFrictMat_NormalInelasticityPhys);
+REGISTER_SERIALIZABLE(Ip2_2xNormalInelasticMat_NormalInelasticityPhys);
=== modified file 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp 2010-04-22 16:40:45 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_CapillaryPhys.cpp 2010-06-04 15:26:30 +0000
@@ -41,13 +41,11 @@
Real fb = sdec2->frictionAngle;
Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);//harmonic average of two stiffnesses
Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);//harmonic average of two stiffnesses with ks=V*kn for each sphere
- contactPhysics->initialKn = Kn;
- contactPhysics->initialKs = Ks;
contactPhysics->frictionAngle = std::min(fa,fb);
contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
contactPhysics->prevNormal = interactionGeometry->normal;
- contactPhysics->kn = contactPhysics->initialKn;
- contactPhysics->ks = contactPhysics->initialKs;
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
}
}
};
=== modified file 'pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 2010-06-03 17:49:44 +0000
+++ pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 2010-06-04 15:26:30 +0000
@@ -8,7 +8,7 @@
#include<yade/pkg-dem/NormalInelasticityLaw.hpp>
-#include<yade/pkg-dem/CohFrictMat.hpp>
+#include<yade/pkg-dem/NormalInelasticMat.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
@@ -21,72 +21,76 @@
const Real& dt = scene->dt;
- int id1 = contact->getId1();
- int id2 = contact->getId2();
-// cout << "contact entre " << id1 << " et " << id2 << " reel ? " << contact->isReal() << endl;
-
- State* de1 = Body::byId(id1,scene)->state.get();
- State* de2 = Body::byId(id2,scene)->state.get();
- ScGeom* currentContactGeometry = YADE_CAST<ScGeom*>(iG.get());
- NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (iP.get());//remplace par :
-
- Vector3r& shearForce = currentContactPhysics->shearForce;
-
- if (contact->isFresh(scene))
- {
- shearForce = Vector3r::Zero();
- currentContactPhysics->previousun=0.0;
- currentContactPhysics->previousFn=0.0;
- currentContactPhysics->unMax=0.0;
- }
-
-
- Real previousun = currentContactPhysics->previousun;
- Real previousFn = currentContactPhysics->previousFn;
- Real kn = currentContactPhysics->kn;
-// cout << "Valeur de kn : " << kn << endl;
- Real Fn; // la valeur de Fn qui va etre calculee selon différentes manieres puis affectee
-
- Real un = currentContactGeometry->penetrationDepth; // compte positivement lorsq vraie interpenetration
-// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << previousun << " et previousFn =" << previousFn << endl;
- if(un > currentContactPhysics->unMax) // on est en charge, et au delà et sur la "droite principale"
- {
- Fn = kn*un;
- currentContactPhysics->unMax = std::abs(un);
-// cout << "je suis dans le calcul normal " << endl;
- }
- else
- {
- Fn = previousFn + coeff_dech * kn * (un-previousun); // Calcul incrémental de la nvlle force
-// cout << "je suis dans l'autre calcul" << endl;
- if(std::abs(Fn) > std::abs(kn * un)) // verif qu'on ne depasse la courbe
- Fn = kn*un;
- if(Fn < 0.0 ) // verif qu'on reste positif FIXME ATTENTION ON S'EST FICHU DU NORMAL ADHESION !!!!
- {Fn = 0;
-// cout << "j'ai corrige pour ne pas etre negatif" << endl;
- }
- }
-
- currentContactPhysics->normalForce = Fn*currentContactGeometry->normal;
-// cout << "Fn appliquee " << Fn << endl << endl;
- // actualisation :
- currentContactPhysics->previousFn = Fn;
- currentContactPhysics->previousun = un;
-
- if ( un < 0 )
- { // BREAK due to tension
-
- //currentContactPhysics->SetBreakingState();
-
-
- scene->interactions->requestErase(contact->getId1(),contact->getId2());
- // probably not useful anymore
+ int id1 = contact->getId1();
+ int id2 = contact->getId2();
+// cout << "contact entre " << id1 << " et " << id2;
+
+ 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());
+ NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (iP.get());
+
+ Vector3r& shearForce = currentContactPhysics->shearForce;
+
+ if (contact->isFresh(scene))
+ {
+ shearForce = Vector3r::Zero();
+ currentContactPhysics->previousun=0.0;
+ currentContactPhysics->previousFn=0.0;
+ currentContactPhysics->unMax=0.0;
+ }
+
+
+ // *** Computation of normal Force : depends of the history *** //
+ if( currentContactGeometry->penetrationDepth < currentContactPhysics->unMax)
+ currentContactPhysics->kn *= Mat1->coeff_dech;
+
+ 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;
+ Real un = currentContactGeometry->penetrationDepth; // compte positivement lorsq vraie interpenetration
+
+// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << currentContactPhysics->previousun << " et previousFn =" << currentContactPhysics->previousFn << endl;
+
+ if(un >= currentContactPhysics->unMax) // on est en charge, et au delà et sur la "droite principale"
+ {
+ Fn = currentContactPhysics->knLower*un;
+ 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
+ {
+ currentContactPhysics->kn = currentContactPhysics->knLower* Mat1->coeff_dech;
+ Fn = currentContactPhysics->previousFn + currentContactPhysics->kn * (un-currentContactPhysics->previousun); // Calcul incrémental de la nvlle force
+// cout << "je suis dans l'autre calcul" << endl;
+ if(std::abs(Fn) > std::abs(currentContactPhysics->knLower * un)) // verif qu'on ne depasse la courbe
+ {
+ Fn = currentContactPhysics->knLower*un;
+// cout << "j'etais dans l'autre calcul mais j'ai corrige pour ne pas depasser la limite" << endl;
+ }
+ if(Fn < 0.0 ) // verif qu'on reste positif
+ {
+ Fn = 0;
+// cout << "j'ai corrige pour ne pas etre negatif" << endl;
+ }
+ }
+ currentContactPhysics->normalForce = Fn*currentContactGeometry->normal;
+// cout << "Fn appliquee " << Fn << endl << endl;
+ // actualisation :
+ currentContactPhysics->previousFn = Fn;
+ currentContactPhysics->previousun = un;
+// *** End of computation of normal force *** //
+
+
+ if ( un < 0 )
+ { // BREAK due to tension
+ scene->interactions->requestErase(contact->getId1(),contact->getId2());
+ // probably not useful anymore
currentContactPhysics->normalForce = Vector3r::Zero();
currentContactPhysics->shearForce = Vector3r::Zero();
- }
- else
- {
-
+ }
+ else
+ {
Vector3r axis;
Real angle;
@@ -138,8 +142,6 @@
if ( Fs > maxFs )
{
- currentContactPhysics->SetBreakingState();
-
maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
maxFs = maxFs / Fs;
=== modified file 'pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp 2010-06-03 17:49:44 +0000
+++ pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp 2010-06-04 15:26:30 +0000
@@ -26,8 +26,7 @@
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity,
LawFunctor,
- "Contact law including cohesion, moment transfer and inelastic compression behaviour\n\n This contact Law is inspired by :yref:`CohesiveFrictionalContactLaw` (inspired itselve directly from the work of Plassiard & Belheine, see the corresponding articles in [DeghmReport2006]_ for example).\n\n It allows so to set moments, cohesion, tension limit and (that's the difference) inelastic unloadings in compression between bodies. All that concerned brokenBodies (this flag and the erosionactivated one) and the useless 'iter' has been suppressed.\n\n The Relationsships corresponding are Ip2_2xCohFrictMat_NormalInelasticityPhys, where the rigidities, the friction angles (with their tan()), and the orientations of the interactions are calculated. No more cohesion and tension limits are computed for all the interactions.\n\n To use it you should also use :\n- :yref:`CohFrictMat` for bodies, with *isCohesive* = 1 (A verifier ce dernier point)\n- :yref:`Ip2_2xCohFrictMat_NormalInelasticityPhys` (=> which involves interactions of :yref:`NormalInelasticityPhys` type).\n\n The effects of this law are illustrated in scripts/NormalInelasticityTest.py",
- ((Real,coeff_dech,1.0,"=kn(unload) / kn(load)"))
+ "Contact law including cohesion, moment transfer and inelastic compression behaviour\n\n For what concerns moment transfer it is inspired indirectly by the work of Plassiard & Belheine, see the corresponding articles in [DeghmReport2006]_ for example).\n\n It allows moreover to set inelastic unloadings in compression between bodies.\n\n The Relationsships corresponding are Ip2_2xNormalInelasticMat_NormalInelasticityPhys, where the rigidities, the friction angles (with their tan()), and the orientations of the interactions are calculated.\n\n To use it you should also use :\n- :yref:`NormalInelasticMat` \n- :yref:`Ip2_2xNormalInelasticMat_NormalInelasticityPhys` (=> which involves interactions of :yref:`NormalInelasticityPhys` type).\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"))
);
=== modified file 'pkg/dem/PreProcessor/SimpleShear.cpp'
--- pkg/dem/PreProcessor/SimpleShear.cpp 2010-06-02 16:40:46 +0000
+++ pkg/dem/PreProcessor/SimpleShear.cpp 2010-06-04 15:26:30 +0000
@@ -15,7 +15,7 @@
#include <yade/pkg-dem/CohFrictMat.hpp>
#include<yade/pkg-dem/NormalInelasticityLaw.hpp>
-#include <yade/pkg-dem/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp>
+#include <yade/pkg-dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp>
#include<yade/pkg-dem/GlobalStiffnessTimeStepper.hpp>
// #include <yade/pkg-dem/PositionRecorder.hpp>
@@ -260,7 +260,7 @@
interactionGeometryDispatcher->add("Ig2_Box_Sphere_ScGeom");
shared_ptr<InteractionPhysicsDispatcher> interactionPhysicsDispatcher(new InteractionPhysicsDispatcher);
- shared_ptr<InteractionPhysicsFunctor> CL1Rel(new Ip2_2xCohFrictMat_NormalInelasticityPhys);
+ shared_ptr<InteractionPhysicsFunctor> CL1Rel(new Ip2_2xNormalInelasticMat_NormalInelasticityPhys);
interactionPhysicsDispatcher->add(CL1Rel);
shared_ptr<BoundDispatcher> boundDispatcher = shared_ptr<BoundDispatcher>(new BoundDispatcher);
=== modified file 'py/system.py'
--- py/system.py 2010-05-02 16:02:29 +0000
+++ py/system.py 2010-06-04 15:26:30 +0000
@@ -141,6 +141,10 @@
'Ip2_SimleViscoelasticMat_SimpleViscoelasticMat_SimpleViscoelasticPhys':'Ip2_ViscElMat_ViscElMat_ViscElPhys', # Fri Apr 9 19:28:48 2010, vaclav@flux
'MomentEngine':'TorqueEngine', # Sun May 2 16:09:34 2010, vaclav@flux
'JumpChangeSe3':'StepDisplacer', # Sun May 2 16:14:21 2010, vaclav@flux
+ 'ContactLaw1':'SomeBetterName', # Fri Jun 4 15:35:38 2010, jduriez@c1solimara-l
+ 'Ip2_2xCohFrictMat_NormalInelasticityPhys':'Ip2_2xNormalInelasticMat_NormalInelasticityPhys', # Fri Jun 4 15:36:41 2010, jduriez@c1solimara-l
+ 'Ip2_2xCohFrictMat_NormalInelasticityPhys':'Ip2_2xNormalInelasticMat_NormalInelasticityPhys', # Fri Jun 4 15:37:01 2010, jduriez@c1solimara-l
+ 'Ip2_2xCohFrictMat_NormalInelasticityPhys':'Ip2_2xNormalInelasticMat_NormalInelasticityPhys', # Fri Jun 4 15:37:16 2010, jduriez@c1solimara-l
### END_RENAMED_CLASSES_LIST ### (do not delete this line; scripts/rename-class.py uses it
}
=== modified file 'scripts/NormalInelasticityTest.py'
--- scripts/NormalInelasticityTest.py 2010-06-03 17:49:44 +0000
+++ scripts/NormalInelasticityTest.py 2010-06-04 15:26:30 +0000
@@ -7,7 +7,7 @@
from yade import plot
#Def of the material which will be used
-O.materials.append(CohFrictMat(density=2600,young=4.0e9,poisson=.04,frictionAngle=.6,label='Materiau1'))
+O.materials.append(NormalInelasticMat(density=2600,young=4.0e9,poisson=.04,frictionAngle=.6,coeff_dech=3.0,label='Materiau1'))
#Def of the bodies of the simulations : 2 spheres, with names which will be useful after
O.bodies.append(utils.sphere([0,0,0], 1, dynamic=False, wire=False, color=None, highlight=False)) #'Materiau1', as the latest material defined will be used
@@ -25,7 +25,7 @@
InteractionDispatchers(
[Ig2_Sphere_Sphere_ScGeom()],
[Ip2_2xCohFrictMat_NormalInelasticityPhys()],
- [Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity(coeff_dech=3)]
+ [Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity()]
),
PeriodicPythonRunner(iterPeriod=1,command='letMove()')
]
