yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11053
[Branch ~yade-pkg/yade/git-trunk] Rev 4077: Move all NormalInelastic files into NormalInelasticPM
------------------------------------------------------------
revno: 4077
committer: Anton Gladky <gladky.anton@xxxxxxxxx>
timestamp: Mon 2014-07-14 22:08:34 +0200
message:
Move all NormalInelastic files into NormalInelasticPM
removed:
pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp
pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp
pkg/dem/NormalInelasticMat.cpp
pkg/dem/NormalInelasticMat.hpp
pkg/dem/NormalInelasticityLaw.cpp
pkg/dem/NormalInelasticityLaw.hpp
pkg/dem/NormalInelasticityPhys.cpp
pkg/dem/NormalInelasticityPhys.hpp
added:
pkg/dem/NormalInelasticPM.cpp
pkg/dem/NormalInelasticPM.hpp
modified:
pkg/dem/SimpleShear.cpp
--
lp:yade
https://code.launchpad.net/~yade-pkg/yade/git-trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-pkg/yade/git-trunk/+edit-subscription
=== removed file 'pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp'
--- pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 2011-08-26 09:49:25 +0000
+++ pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,69 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* 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"Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp"
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/NormalInelasticityPhys.hpp>
-#include<yade/pkg/dem/NormalInelasticMat.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-
-
-
-void Ip2_2xNormalInelasticMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // NormalInelasticMat
- , const shared_ptr<Material>& b2 // NormalInelasticMat
- , const shared_ptr<Interaction>& interaction)
-{
- NormalInelasticMat* sdec1 = static_cast<NormalInelasticMat*>(b1.get());
- NormalInelasticMat* sdec2 = static_cast<NormalInelasticMat*>(b2.get());
- ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
-
-
- if(geom) // so it is ScGeom - NON PERMANENT LINK
- {
- if(!interaction->phys)
- {
-//std::cerr << " isNew, id1: " << interaction->getId1() << " id2: " << interaction->getId2() << "\n";
- interaction->phys = shared_ptr<NormalInelasticityPhys>(new NormalInelasticityPhys());
- NormalInelasticityPhys* contactPhysics = YADE_CAST<NormalInelasticityPhys*>(interaction->phys.get());
-
- Real Ea = sdec1->young;
- Real Eb = sdec2->young;
- Real Va = sdec1->poisson;
- Real Vb = sdec2->poisson;
- Real Ra = geom->radius1;
- Real Rb = geom->radius2;
- Real fa = sdec1->frictionAngle;
- Real fb = sdec2->frictionAngle;
-
- Real Kn = 2.0*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);//harmonic average of two stiffnesses
-
- Real Ks = 2.0*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
-
- // Jean-Patrick Plassiard, Noura Belheine, Frederic Victor Donze, "A Spherical Discrete Element Model: calibration procedure and incremental response", DOI: 10.1007/s10035-009-0130-x
-
- Real Kr = betaR*std::pow((Ra+Rb)/2.0,2)*Ks;
-
- contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
- contactPhysics->forMaxMoment = 1.0*(Ra+Rb)/2.0; // 1.0 corresponding to ethaR which I don't know exactly where to define as a parameter...
-
- // Lot of suppress here around (>) r2276.
- contactPhysics->knLower = Kn;
- contactPhysics->kn = Kn;
- contactPhysics->ks = Ks;
- contactPhysics->kr = Kr;
- }
-
- }
-
-};
-YADE_PLUGIN((Ip2_2xNormalInelasticMat_NormalInelasticityPhys));
-
-
-
=== removed file 'pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp'
--- pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp 2011-05-04 10:19:03 +0000
+++ pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,38 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* 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/Dispatching.hpp>
-#include<yade/pkg/dem/NormalInelasticMat.hpp>
-
-/*! \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"
- */
-
-class Ip2_2xNormalInelasticMat_NormalInelasticityPhys : public IPhysFunctor
-{
- public :
-
- virtual void go( const shared_ptr<Material>& b1,
- const shared_ptr<Material>& b2,
- const shared_ptr<Interaction>& interaction);
-
- FUNCTOR2D(NormalInelasticMat,NormalInelasticMat);
- YADE_CLASS_BASE_DOC_ATTRS(Ip2_2xNormalInelasticMat_NormalInelasticityPhys,
- IPhysFunctor,
- "The RelationShips for using :yref:`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"))
- );
-};
-
-REGISTER_SERIALIZABLE(Ip2_2xNormalInelasticMat_NormalInelasticityPhys);
-
-
=== removed file 'pkg/dem/NormalInelasticMat.cpp'
--- pkg/dem/NormalInelasticMat.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NormalInelasticMat.cpp 1970-01-01 00:00:00 +0000
@@ -1,12 +0,0 @@
-
-#include "NormalInelasticMat.hpp"
-
-
-NormalInelasticMat::~NormalInelasticMat()
-{
-}
-
-YADE_PLUGIN((NormalInelasticMat));
-
-
-
=== removed file 'pkg/dem/NormalInelasticMat.hpp'
--- pkg/dem/NormalInelasticMat.hpp 2011-05-04 12:37:24 +0000
+++ pkg/dem/NormalInelasticMat.hpp 1970-01-01 00:00:00 +0000
@@ -1,32 +0,0 @@
-/*************************************************************************
-* 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 :yref:`coeff_dech<NormalInelasticMat::coeff_dech>`).",
- ((Real,coeff_dech,1.0,,"=kn(unload) / kn(load)"))
- ,
- createIndex();
- );
-/// Indexable
- REGISTER_CLASS_INDEX(NormalInelasticMat,FrictMat);
-};
-
-REGISTER_SERIALIZABLE(NormalInelasticMat);
-
-
=== added file 'pkg/dem/NormalInelasticPM.cpp'
--- pkg/dem/NormalInelasticPM.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/NormalInelasticPM.cpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,198 @@
+/*************************************************************************
+* Copyright (C) 2008 by Jérôme DURIEZ *
+* duriez@xxxxxxxxxxxxxxx *
+* *
+* 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 "NormalInelasticPM.hpp"
+#include <yade/pkg/dem/ScGeom.hpp>
+#include <yade/core/Omega.hpp>
+#include <yade/core/Scene.hpp>
+
+YADE_PLUGIN((NormalInelasticMat)(NormalInelasticityPhys)(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity)(Ip2_2xNormalInelasticMat_NormalInelasticityPhys));
+
+
+void Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
+{
+ int id1 = contact->getId1();
+ int id2 = contact->getId2();
+// cout << "contact entre " << id1 << " et " << id2;
+
+ NormalInelasticMat* Mat1 = static_cast<NormalInelasticMat*>(Body::byId(id1,scene)->material.get());
+ ScGeom6D* geom = YADE_CAST<ScGeom6D*>(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;
+ }
+
+ un = geom->penetrationDepth; // >0 for real penetration
+
+// Check if there is a real overlap or not. The Ig2... seems to let exist interactions with negative un (= no overlap). Such interactions seem then to have to be deleted here.
+ if ( un < 0 )
+ {
+ scene->interactions->requestErase(contact);// this, among other things, resets the interaction : geometry and physics variables (as forces, ...) are reset to defaut values
+ return;
+ }
+
+
+
+// ******** Computation of normal Force : depends of the history ******* //
+
+// cout << " Dans Law2 valeur de kn : " << currentContactPhysics->kn << endl;
+// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << currentContactPhysics->previousun << " et previousFn =" << currentContactPhysics->previousFn << endl;
+
+ if(un >= currentContactPhysics->unMax) // case of virgin load : on the "principal line" (limit state of the (un,Fn) space)
+ {
+ 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 job
+ {
+ currentContactPhysics->kn = currentContactPhysics->knLower* Mat1->coeff_dech;
+ // 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 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;
+ }
+ if(Fn < 0.0 ) // check to stay >=0
+ {
+ Fn = 0;
+// cout << "j'ai corrige pour ne pas etre negatif" << endl;
+ }
+ }
+ currentContactPhysics->normalForce = Fn*geom->normal;
+// cout << "Fn appliquee " << Fn << endl << endl;
+
+ // actualisation :
+ currentContactPhysics->previousFn = Fn;
+ currentContactPhysics->previousun = un;
+
+// *** End of computation of normal force *** //
+
+
+
+// ******** Tangential force ******* //
+ if ( un < 0 )
+ { // BREAK due to tension
+ scene->interactions->requestErase(contact); return;
+ // probably not useful anymore
+ currentContactPhysics->normalForce = Vector3r::Zero();
+ currentContactPhysics->shearForce = Vector3r::Zero();
+ }
+ else
+ {
+ // Correction of previous shear force to take into account the change in normal:
+ shearForce = geom->rotate(currentContactPhysics->shearForce);
+ // Update of shear force corresponding to shear displacement increment:
+ shearForce -= currentContactPhysics->ks * geom->shearIncrement();
+
+ Fs = currentContactPhysics->shearForce.norm();
+ maxFs = std::max((Real) 0,Fn*currentContactPhysics->tangensOfFrictionAngle);
+
+ if ( Fs > maxFs )
+ {
+ maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
+
+ maxFs = maxFs / Fs;
+ if (maxFs>1)
+ cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!" << endl;
+ shearForce *= maxFs;
+ if (Fn<0) currentContactPhysics->normalForce = Vector3r::Zero();
+ }
+
+
+ f = currentContactPhysics->normalForce + shearForce;
+ scene->forces.addForce (id1,-f);
+ scene->forces.addForce (id2, f);
+ scene->forces.addTorque(id1,-(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(f));
+ scene->forces.addTorque(id2, -(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(f));
+// *** End of computation of tangential force *** //
+
+
+// ******** Moment law ******* //
+
+ if(momentRotationLaw)
+ {
+ currentContactPhysics->moment_twist = (geom->getTwist()*currentContactPhysics->kr)*geom->normal ;
+ currentContactPhysics->moment_bending = geom->getBending() * currentContactPhysics->kr;
+
+ moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
+
+// Limitation by plastic threshold of this part of the moment caused by relative twist and bending
+ if (!momentAlwaysElastic)
+ {
+ Real normeMomentMax = currentContactPhysics->forMaxMoment * std::abs(Fn);
+ if(moment.norm()>normeMomentMax)
+ {
+ moment *= normeMomentMax/moment.norm();
+ }
+ }
+ scene->forces.addTorque(id1,-moment);
+ scene->forces.addTorque(id2, moment);
+ }
+// ******** Moment law END ******* //
+ }
+}
+
+
+void Ip2_2xNormalInelasticMat_NormalInelasticityPhys::go( const shared_ptr<Material>& b1 // NormalInelasticMat
+ , const shared_ptr<Material>& b2 // NormalInelasticMat
+ , const shared_ptr<Interaction>& interaction)
+{
+ NormalInelasticMat* sdec1 = static_cast<NormalInelasticMat*>(b1.get());
+ NormalInelasticMat* sdec2 = static_cast<NormalInelasticMat*>(b2.get());
+ ScGeom* geom = YADE_CAST<ScGeom*>(interaction->geom.get());
+
+
+ if(geom) // so it is ScGeom - NON PERMANENT LINK
+ {
+ if(!interaction->phys)
+ {
+//std::cerr << " isNew, id1: " << interaction->getId1() << " id2: " << interaction->getId2() << "\n";
+ interaction->phys = shared_ptr<NormalInelasticityPhys>(new NormalInelasticityPhys());
+ NormalInelasticityPhys* contactPhysics = YADE_CAST<NormalInelasticityPhys*>(interaction->phys.get());
+
+ Real Ea = sdec1->young;
+ Real Eb = sdec2->young;
+ Real Va = sdec1->poisson;
+ Real Vb = sdec2->poisson;
+ Real Ra = geom->radius1;
+ Real Rb = geom->radius2;
+ Real fa = sdec1->frictionAngle;
+ Real fb = sdec2->frictionAngle;
+
+ Real Kn = 2.0*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);//harmonic average of two stiffnesses
+
+ Real Ks = 2.0*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Vb);//harmonic average of two stiffnesses with ks=V*kn for each sphere
+
+ // Jean-Patrick Plassiard, Noura Belheine, Frederic Victor Donze, "A Spherical Discrete Element Model: calibration procedure and incremental response", DOI: 10.1007/s10035-009-0130-x
+
+ Real Kr = betaR*std::pow((Ra+Rb)/2.0,2)*Ks;
+
+ contactPhysics->tangensOfFrictionAngle = std::tan(std::min(fa,fb));
+ contactPhysics->forMaxMoment = 1.0*(Ra+Rb)/2.0; // 1.0 corresponding to ethaR which I don't know exactly where to define as a parameter...
+
+ // Lot of suppress here around (>) r2276.
+ contactPhysics->knLower = Kn;
+ contactPhysics->kn = Kn;
+ contactPhysics->ks = Ks;
+ contactPhysics->kr = Kr;
+ }
+
+ }
+
+};
+
=== added file 'pkg/dem/NormalInelasticPM.hpp'
--- pkg/dem/NormalInelasticPM.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/NormalInelasticPM.hpp 2014-07-14 20:08:34 +0000
@@ -0,0 +1,114 @@
+/*************************************************************************
+* 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>
+#include <yade/pkg/dem/FrictPhys.hpp>
+#include <yade/pkg/common/Dispatching.hpp>
+#include <yade/pkg/dem/ScGeom.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 :yref:`coeff_dech<NormalInelasticMat::coeff_dech>`).",
+ ((Real,coeff_dech,1.0,,"=kn(unload) / kn(load)"))
+ ,
+ createIndex();
+ );
+/// Indexable
+ REGISTER_CLASS_INDEX(NormalInelasticMat,FrictMat);
+};
+
+REGISTER_SERIALIZABLE(NormalInelasticMat);
+
+
+class NormalInelasticityPhys : public FrictPhys
+{
+ public :
+ virtual ~NormalInelasticityPhys() {};
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
+ "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"))
+ ((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 interaction"))
+ ((Real,knLower,0.0,,"the stifness corresponding to a virgin load for example"))
+ // internal attributes
+ ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
+ ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
+ ,
+ createIndex();
+ );
+ REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhys);
+};
+
+REGISTER_SERIALIZABLE(NormalInelasticityPhys);
+
+
+class Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
+{
+ private :
+ Vector3r moment // the part of the contact torque of the interaction due to relative rotations (a first part is due to contact forces)
+ ,f// contact force
+ ;
+ Real Fn // value of normal force in the interaction
+ ,Fs // shear force
+ ,maxFs; // maximum value of shear force according to Coulomb-like criterion
+ Real un; // value of interpenetration in the interaction
+ public :
+ virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+
+ FUNCTOR2D(ScGeom,NormalInelasticityPhys);
+
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity,
+ LawFunctor,
+ "Contact law used to simulate granular filler in rock joints [Duriez2009a]_, [Duriez2011]_. 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 :yref:`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 examples/normalInelasticity-test.py",
+ ((bool,momentRotationLaw,true,,"boolean, true=> computation of a torque (against relative rotation) exchanged between particles"))
+ ((bool,momentAlwaysElastic,false,,"boolean, true=> the part of the contact torque (caused by relative rotations, which is computed only if momentRotationLaw..) is not limited by a plastic threshold"))
+ ,
+ moment=Vector3r::Zero();
+ f=Vector3r::Zero();
+ Fn=0.0;
+ Fs=0.0;
+ maxFs=0.0;
+ un=0.0;
+ );
+
+};
+
+REGISTER_SERIALIZABLE(Law2_ScGeom6D_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"
+ */
+
+class Ip2_2xNormalInelasticMat_NormalInelasticityPhys : public IPhysFunctor
+{
+ public :
+
+ virtual void go( const shared_ptr<Material>& b1,
+ const shared_ptr<Material>& b2,
+ const shared_ptr<Interaction>& interaction);
+
+ FUNCTOR2D(NormalInelasticMat,NormalInelasticMat);
+ YADE_CLASS_BASE_DOC_ATTRS(Ip2_2xNormalInelasticMat_NormalInelasticityPhys,
+ IPhysFunctor,
+ "The RelationShips for using :yref:`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"))
+ );
+};
+
+REGISTER_SERIALIZABLE(Ip2_2xNormalInelasticMat_NormalInelasticityPhys);
+
=== removed file 'pkg/dem/NormalInelasticityLaw.cpp'
--- pkg/dem/NormalInelasticityLaw.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/NormalInelasticityLaw.cpp 1970-01-01 00:00:00 +0000
@@ -1,153 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* 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<yade/pkg/dem/NormalInelasticityLaw.hpp>
-
-#include<yade/pkg/dem/NormalInelasticMat.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-YADE_PLUGIN((Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity));
-
-
-
-void Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
-{
- int id1 = contact->getId1();
- int id2 = contact->getId2();
-// cout << "contact entre " << id1 << " et " << id2;
-
- NormalInelasticMat* Mat1 = static_cast<NormalInelasticMat*>(Body::byId(id1,scene)->material.get());
- ScGeom6D* geom = YADE_CAST<ScGeom6D*>(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;
- }
-
- un = geom->penetrationDepth; // >0 for real penetration
-
-// Check if there is a real overlap or not. The Ig2... seems to let exist interactions with negative un (= no overlap). Such interactions seem then to have to be deleted here.
- if ( un < 0 )
- {
- scene->interactions->requestErase(contact);// this, among other things, resets the interaction : geometry and physics variables (as forces, ...) are reset to defaut values
- return;
- }
-
-
-
-// ******** Computation of normal Force : depends of the history ******* //
-
-// cout << " Dans Law2 valeur de kn : " << currentContactPhysics->kn << endl;
-// cout << "un = " << un << " alors que unMax = "<< currentContactPhysics->unMax << " et previousun = " << currentContactPhysics->previousun << " et previousFn =" << currentContactPhysics->previousFn << endl;
-
- if(un >= currentContactPhysics->unMax) // case of virgin load : on the "principal line" (limit state of the (un,Fn) space)
- {
- 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 job
- {
- currentContactPhysics->kn = currentContactPhysics->knLower* Mat1->coeff_dech;
- // 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 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;
- }
- if(Fn < 0.0 ) // check to stay >=0
- {
- Fn = 0;
-// cout << "j'ai corrige pour ne pas etre negatif" << endl;
- }
- }
- currentContactPhysics->normalForce = Fn*geom->normal;
-// cout << "Fn appliquee " << Fn << endl << endl;
-
- // actualisation :
- currentContactPhysics->previousFn = Fn;
- currentContactPhysics->previousun = un;
-
-// *** End of computation of normal force *** //
-
-
-
-// ******** Tangential force ******* //
- if ( un < 0 )
- { // BREAK due to tension
- scene->interactions->requestErase(contact); return;
- // probably not useful anymore
- currentContactPhysics->normalForce = Vector3r::Zero();
- currentContactPhysics->shearForce = Vector3r::Zero();
- }
- else
- {
- // Correction of previous shear force to take into account the change in normal:
- shearForce = geom->rotate(currentContactPhysics->shearForce);
- // Update of shear force corresponding to shear displacement increment:
- shearForce -= currentContactPhysics->ks * geom->shearIncrement();
-
- Fs = currentContactPhysics->shearForce.norm();
- maxFs = std::max((Real) 0,Fn*currentContactPhysics->tangensOfFrictionAngle);
-
- if ( Fs > maxFs )
- {
- maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
-
- maxFs = maxFs / Fs;
- if (maxFs>1)
- cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!" << endl;
- shearForce *= maxFs;
- if (Fn<0) currentContactPhysics->normalForce = Vector3r::Zero();
- }
-
-
- f = currentContactPhysics->normalForce + shearForce;
- scene->forces.addForce (id1,-f);
- scene->forces.addForce (id2, f);
- scene->forces.addTorque(id1,-(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(f));
- scene->forces.addTorque(id2, -(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(f));
-// *** End of computation of tangential force *** //
-
-
-// ******** Moment law ******* //
-
- if(momentRotationLaw)
- {
- currentContactPhysics->moment_twist = (geom->getTwist()*currentContactPhysics->kr)*geom->normal ;
- currentContactPhysics->moment_bending = geom->getBending() * currentContactPhysics->kr;
-
- moment = currentContactPhysics->moment_twist + currentContactPhysics->moment_bending;
-
-// Limitation by plastic threshold of this part of the moment caused by relative twist and bending
- if (!momentAlwaysElastic)
- {
- Real normeMomentMax = currentContactPhysics->forMaxMoment * std::abs(Fn);
- if(moment.norm()>normeMomentMax)
- {
- moment *= normeMomentMax/moment.norm();
- }
- }
- scene->forces.addTorque(id1,-moment);
- scene->forces.addTorque(id2, moment);
- }
-// ******** Moment law END ******* //
- }
-}
-
-
-
=== removed file 'pkg/dem/NormalInelasticityLaw.hpp'
--- pkg/dem/NormalInelasticityLaw.hpp 2014-05-28 08:17:25 +0000
+++ pkg/dem/NormalInelasticityLaw.hpp 1970-01-01 00:00:00 +0000
@@ -1,51 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* 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/Dispatching.hpp>
-#include<yade/pkg/dem/ScGeom.hpp>
-#include<yade/pkg/dem/NormalInelasticityPhys.hpp>
-#include <set>
-#include <boost/tuple/tuple.hpp>
-
-
-
-class Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
-{
- private :
- Vector3r moment // the part of the contact torque of the interaction due to relative rotations (a first part is due to contact forces)
- ,f// contact force
- ;
- Real Fn // value of normal force in the interaction
- ,Fs // shear force
- ,maxFs; // maximum value of shear force according to Coulomb-like criterion
- Real un; // value of interpenetration in the interaction
- public :
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
-
- FUNCTOR2D(ScGeom,NormalInelasticityPhys);
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity,
- LawFunctor,
- "Contact law used to simulate granular filler in rock joints [Duriez2009a]_, [Duriez2011]_. 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 :yref:`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 examples/normalInelasticity-test.py",
- ((bool,momentRotationLaw,true,,"boolean, true=> computation of a torque (against relative rotation) exchanged between particles"))
- ((bool,momentAlwaysElastic,false,,"boolean, true=> the part of the contact torque (caused by relative rotations, which is computed only if momentRotationLaw..) is not limited by a plastic threshold"))
- ,
- moment=Vector3r::Zero();
- f=Vector3r::Zero();
- Fn=0.0;
- Fs=0.0;
- maxFs=0.0;
- un=0.0;
- );
-
-};
-
-REGISTER_SERIALIZABLE(Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity);
-
=== removed file 'pkg/dem/NormalInelasticityPhys.cpp'
--- pkg/dem/NormalInelasticityPhys.cpp 2010-10-13 16:23:08 +0000
+++ pkg/dem/NormalInelasticityPhys.cpp 1970-01-01 00:00:00 +0000
@@ -1,18 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* 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 "NormalInelasticityPhys.hpp"
-
-
-NormalInelasticityPhys::~NormalInelasticityPhys()
-{
-}
-
-
-YADE_PLUGIN((NormalInelasticityPhys));
-
=== removed file 'pkg/dem/NormalInelasticityPhys.hpp'
--- pkg/dem/NormalInelasticityPhys.hpp 2011-01-24 14:45:35 +0000
+++ pkg/dem/NormalInelasticityPhys.hpp 1970-01-01 00:00:00 +0000
@@ -1,38 +0,0 @@
-/*************************************************************************
-* Copyright (C) 2008 by Jérôme DURIEZ *
-* duriez@xxxxxxxxxxxxxxx *
-* *
-* 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/dem/FrictPhys.hpp>
-
-
-class NormalInelasticityPhys : public FrictPhys
-{
- public :
- virtual ~NormalInelasticityPhys();
-
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(NormalInelasticityPhys,FrictPhys,
- "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"))
- ((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 interaction"))
- ((Real,knLower,0.0,,"the stifness corresponding to a virgin load for example"))
- // internal attributes
- ((Vector3r,moment_twist,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Twist moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
- ((Vector3r,moment_bending,Vector3r(0,0,0),(Attr::noSave | Attr::readonly),"Bending moment. Defined here, being initialized as it should be, to be used in :yref:`Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity`"))
- ,
- createIndex();
- );
- REGISTER_CLASS_INDEX(NormalInelasticityPhys,FrictPhys);
-};
-
-REGISTER_SERIALIZABLE(NormalInelasticityPhys);
-
-
=== modified file 'pkg/dem/SimpleShear.cpp'
--- pkg/dem/SimpleShear.cpp 2014-07-02 16:18:24 +0000
+++ pkg/dem/SimpleShear.cpp 2014-07-14 20:08:34 +0000
@@ -11,9 +11,7 @@
#include "SimpleShear.hpp"
-#include <yade/pkg/dem/NormalInelasticMat.hpp>
-#include<yade/pkg/dem/NormalInelasticityLaw.hpp>
-#include <yade/pkg/dem/Ip2_2xNormalInelasticMat_NormalInelasticityPhys.hpp>
+#include <yade/pkg/dem/NormalInelasticPM.hpp>
#include<yade/pkg/dem/GlobalStiffnessTimeStepper.hpp>
#include<yade/pkg/common/Aabb.hpp>
