yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #03883
[Branch ~yade-dev/yade/trunk] Rev 2132: - the code concerning the inelastic normal law is now in files named NormalInelasticityLaw.* inst...
------------------------------------------------------------
revno: 2132
committer: jduriez <jduriez@c1solimara-l>
branch nick: trunk
timestamp: Thu 2010-04-08 11:14:11 +0200
message:
- the code concerning the inelastic normal law is now in files named NormalInelasticityLaw.* instead of Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.*
The class itself has still the same name (Law2...), but I changed the name of the files for practical reasons, as it seem other people do (files HertzMindlin,
CohesiveFrictionalPM for ex). Tell me if this is problematic
- Changes in "include" of related files
- Very beginning of changing this law in a functor instead of an Engine (draft lines which are yet commented)
- File Ip2_FrictMat_FrictMat_FrictPhys.cpp appear to have been modified for "historical" reasons on my side, but it is not the case
removed:
pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.cpp
pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp
added:
pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp
pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp
modified:
pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp
pkg/dem/Engine/PartialEngine/KinemCNLEngine.hpp
pkg/dem/Engine/PartialEngine/KinemCNSEngine.hpp
pkg/dem/PreProcessor/SimpleShear.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 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-03-27 22:18:10 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-04-08 09:14:11 +0000
@@ -34,13 +34,13 @@
Real Vb = mat2->poisson;
Real Da,Db; Vector3r normal;
- //FIXME : dynamic casts here???!!!
+ //FIXME : dynamic casts here???!!!
ScGeom* scg=dynamic_cast<ScGeom*>(interaction->interactionGeometry.get());
Dem3DofGeom* d3dg=dynamic_cast<Dem3DofGeom*>(interaction->interactionGeometry.get());
if(scg){Da=scg->radius1; Db=scg->radius2; normal=scg->normal;}
else if(d3dg){Da=d3dg->refR1>0?d3dg->refR1:2*d3dg->refR2; Db=d3dg->refR2>0?d3dg->refR2:d3dg->refR1; normal=d3dg->normal;}
else throw runtime_error("Ip2_FrictMat_FrictMat_FrictPhys: geometry is neither ScGeom nor Dem3DofGeom");
- //harmonic average of the two stiffnesses when (Di.Ei/2) is the stiffness of sphere "i"
+ //harmonic average of the two stiffnesses when (Di.Ei/2) is the stiffness of sphere "i"
Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);
//same for shear stiffness
Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);
=== removed file 'pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.cpp'
--- pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.cpp 2010-03-23 14:19:56 +0000
+++ pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.cpp 1970-01-01 00:00:00 +0000
@@ -1,295 +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/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp>
-
-#include<yade/pkg-dem/CohesiveFrictionalMat.hpp>
-#include<yade/pkg-dem/ScGeom.hpp>
-#include<yade/pkg-dem/NormalInelasticityPhys.hpp>
-#include<yade/core/Omega.hpp>
-#include<yade/core/Scene.hpp>
-
-YADE_PLUGIN((Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity));
-
-
-Vector3r translation_vect (0.10,0,0);
-
-
-void Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity::action()
-{
- shared_ptr<BodyContainer>& bodies = scene->bodies;
-
- Real dt = Omega::instance().getTimeStep();
-
- InteractionContainer::iterator ii = scene->interactions->begin();
- InteractionContainer::iterator iiEnd = scene->interactions->end();
- int nbreInteracTot=0;
- int nbreInteracMomPlastif=0;
- for ( ; ii!=iiEnd ; ++ii )
- {
- //if ((*ii)->interactionGeometry && (*ii)->interactionPhysics)
- if ((*ii)->isReal())
- {
- nbreInteracTot++;
- const shared_ptr<Interaction>& contact = *ii;
- int id1 = contact->getId1();
- int id2 = contact->getId2();
-// cout << "contact entre " << id1 << " et " << id2 << " reel ? " << contact->isReal() << endl;
- if ( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) )
- continue; // skip other groups,
-
-// CohesiveFrictionalMat* de1 = YADE_CAST<CohesiveFrictionalMat*>((*bodies)[id1]->physicalParameters.get());
- State* de1 = Body::byId(id1,scene)->state.get();
- State* de2 = Body::byId(id2,scene)->state.get();
- ScGeom* currentContactGeometry = YADE_CAST<ScGeom*>(contact->interactionGeometry.get());
- NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (contact->interactionPhysics.get());
-
- 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;
- 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
- currentContactPhysics->normalForce = Vector3r::ZERO;
- currentContactPhysics->shearForce = Vector3r::ZERO;
-
- //return;
- // else
- // currentContactPhysics->normalForce = -currentContactPhysics->normalAdhesion*currentContactGeometry->normal;
- }
- else
- {
-
- Vector3r axis;
- Real angle;
-
- /// Here is the code with approximated rotations ///
-
- axis = currentContactPhysics->prevNormal.Cross(currentContactGeometry->normal);
- shearForce -= shearForce.Cross(axis);
- angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angVel+de2->angVel);
- axis = angle*currentContactGeometry->normal;
- shearForce -= shearForce.Cross(axis);
-
- /// Here is the code with exact rotations ///
-
- // Quaternionr q;
- //
- // axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
- // angle = acos(currentContactGeometry->normal.dot(currentContactPhysics->prevNormal));
- // q.fromAngleAxis(angle,axis);
- //
- // currentContactPhysics->shearForce = currentContactPhysics->shearForce*q;
- //
- // angle = dt*0.5*currentContactGeometry->normal.dot(de1->angVel+de2->angVel);
- // axis = currentContactGeometry->normal;
- // 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)
- 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_));
- Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.Dot(relativeVelocity)*currentContactGeometry->normal;
- Vector3r shearDisplacement = shearVelocity*dt;
-
-
-///////////////////////// CREEP START (commented out) ///////////
-// Real viscosity = 300000.0;
-// shearForce -= currentContactPhysics->ks*(shearDisplacement + shearForce*dt/viscosity);
-
- shearForce -= currentContactPhysics->ks*shearDisplacement;
-///////////////////////// CREEP END /////////////////////////////
-
- // cerr << "shearForce = " << shearForce << endl;
- Real maxFs = 0;
- Real Fs = currentContactPhysics->shearForce.Length();
- maxFs = std::max((Real) 0,Fn*currentContactPhysics->tangensOfFrictionAngle);
-
- if ( Fs > maxFs )
- {
-
-
- currentContactPhysics->SetBreakingState();
-
- // maxFs = currentContactPhysics->shearAdhesion;
- // if (!currentContactPhysics->cohesionDisablesFriction && un>0)
- // maxFs += Fn * currentContactPhysics->tangensOfFrictionAngle;
-
- maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
-
- //cerr << "currentContactPhysics->tangensOfFrictionAngle = " << currentContactPhysics->tangensOfFrictionAngle << endl;
- // cerr << "maxFs = " << maxFs << endl;
-
- maxFs = maxFs / Fs;
- // cerr << "maxFs = " << maxFs << endl;
- if (maxFs>1)
- cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
- shearForce *= maxFs;
- if (Fn<0) currentContactPhysics->normalForce = Vector3r::ZERO;
- }
-
-
-
- ////////// PFC3d SlipModel
-
- Vector3r f = currentContactPhysics->normalForce + shearForce;
- // cerr << "currentContactPhysics->normalForce= " << currentContactPhysics->normalForce << endl;
- // cerr << "shearForce " << shearForce << endl;
- // cerr << "f= " << f << endl;
- // it will be some macro( body->physicalActions, ActionType , bodyId )
- scene->forces.addForce (id1,-f);
- scene->forces.addForce (id2,+f);
- scene->forces.addTorque(id1,-c1x.Cross(f));
- scene->forces.addTorque(id2, c2x.Cross(f));
-
-///// /// Moment law ///
-///// if(momentRotationLaw /*&& currentContactPhysics->cohesionBroken == false*/ )
-///// {
-///// Quaternionr delta = ...
-/////
-///// //This indentation is a rewrite of original equations (the two commented lines), should work exactly the same.
-///////Real elasticMoment = currentContactPhysics->kr * std::abs(angle); // positive value (*)
-/////
-///// Real angle_twist(angle * axis.Dot(currentContactGeometry->normal) );
-///// Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
-//////* no creep
-/////*/
-///// Vector3r moment_twist(axis_twist * currentContactPhysics->kr);
-/////
-///// Vector3r axis_bending(angle*axis - axis_twist);
-///// Vector3r moment_bending(axis_bending * currentContactPhysics->kr);
-/////
-///// /*
-///// // We cannot have insanely big moment, so we take a min value of ELASTIC and PLASTIC moment.
-///// Real avgRadius = 0.5 * (currentContactGeometry->radius1 + currentContactGeometry->radius2);
-///// // FIXME - elasticRollingLimit is currently assumed to be 1.0
-///// Real plasticMoment = currentContactPhysics->elasticRollingLimit * avgRadius * std::abs(Fn); // positive value (*)
-///// Real moment(std::min(elasticMoment, plasticMoment)); // RESULT
-///// */
-/////
-///////Vector3r moment = axis * elasticMoment * (angle<0.0?-1.0:1.0); // restore sign. (*)
-/////
-///// Vector3r moment = moment_twist + moment_bending;
-/////
-///// static_cast<Momentum*>( scene->physicalActions->find( id1 , actionMomentum->getClassIndex() ).get() )->momentum += moment;
-///// static_cast<Momentum*>( scene->physicalActions->find( id2 , actionMomentum->getClassIndex() ).get() )->momentum -= moment;
-///// }
-///// /// Moment law END ///
-
- /// Moment law ///
- if(momentRotationLaw)
- {
- {// updates only orientation of contact (local coordinate system)
- Vector3r axis = currentContactPhysics->prevNormal.UnitCross(currentContactGeometry->normal);
- Real angle = unitVectorsAngle(currentContactPhysics->prevNormal,currentContactGeometry->normal);
- Quaternionr align(axis,angle);
- currentContactPhysics->currentContactOrientation = align * currentContactPhysics->currentContactOrientation;
- }
-
- Quaternionr delta( de1->se3.orientation * currentContactPhysics->initialOrientation1.Conjugate() *
- currentContactPhysics->initialOrientation2 * de2->se3.orientation.Conjugate());
-
- Vector3r axis; // axis of rotation - this is the Moment direction UNIT vector.
- Real angle; // angle represents the power of resistant ELASTIC moment
- delta.ToAxisAngle(axis,angle);
- if(angle > Mathr::PI) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
-
- //This indentation is a rewrite of original equations (the two commented lines), should work exactly the same.
-//Real elasticMoment = currentContactPhysics->kr * std::abs(angle); // positive value (*)
-
- Real angle_twist(angle * axis.Dot(currentContactGeometry->normal) );
- Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
- Vector3r moment_twist(axis_twist * currentContactPhysics->kr);
-
- Vector3r axis_bending(angle*axis - axis_twist);
- Vector3r moment_bending(axis_bending * currentContactPhysics->kr);
-
-//Vector3r moment = axis * elasticMoment * (angle<0.0?-1.0:1.0); // restore sign. (*)
-
- Vector3r moment = moment_twist + moment_bending;
-
-// Limitation par seuil plastique
- if (!momentAlwaysElastic)
- {
- Real normeMomentMax = currentContactPhysics->forMaxMoment * std::fabs(Fn);
- Real normeMoment = moment.Length();
- if(normeMoment>normeMomentMax)
- {
- moment *= normeMomentMax/normeMoment;
- nbreInteracMomPlastif++;
- }
- }
- scene->forces.addTorque(id1,-moment);
- scene->forces.addTorque(id2,+moment);
- }
- /// Moment law END ///
-
- currentContactPhysics->prevNormal = currentContactGeometry->normal;
- }
- }
- }
-// cout << " En tout :"<< nbreInteracTot << "interactions (reelles)" << endl;
- //cerr << "ncount= " << ncount << endl;//REMOVE
-// cout << momentAlwaysElastic << endl;
-// cout << "Sur " << nbreInteracTot << " interactions (reelles) " << nbreInteracMomPlastif << " se sont vues corriger leur moment" << endl;
-
-}
-
=== removed file 'pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp'
--- pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp 2010-04-04 22:05:37 +0000
+++ pkg/dem/Engine/GlobalEngine/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp 1970-01-01 00:00:00 +0000
@@ -1,39 +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/LawFunctor.hpp>
-#include<yade/core/GlobalEngine.hpp>
-#include <set>
-#include <boost/tuple/tuple.hpp>
-
-
-
-class Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity : public GlobalEngine
-{
- public :
-// <<<<<<< TREE
- void action();
- YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity,
-// =======
-// void action();
-// YADE_CLASS_BASE_DOC_ATTRS(NormalInelasticityLaw,
-// >>>>>>> MERGE-SOURCE
- GlobalEngine,
- "Contact law including cohesion, moment transfer and inelastic compression behaviour\n\n This contact Law is inspired by :yref:`CohesiveFrictionalContactLaw` (inspired itself directly by the work of Plassiard & Belheine, see 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:`CohesiveFrictionalMat` for bodies, with *isCohesive* = 1 (A verifier ce dernier point)\n- :yref:`Ip2_2xCohFrictMat_NormalInelasticityPhys` (=> which involves interactions of :yref:`NormalInelasticityPhys` type)",
- ((int,sdecGroupMask,1,"?"))
- ((Real,coeff_dech,1.0,"=kn(unload) / kn(load)"))
- ((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);
-
-
=== added file 'pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.cpp 2010-04-08 09:14:11 +0000
@@ -0,0 +1,252 @@
+/*************************************************************************
+* 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/CohesiveFrictionalMat.hpp>
+#include<yade/pkg-dem/ScGeom.hpp>
+#include<yade/pkg-dem/NormalInelasticityPhys.hpp>
+#include<yade/core/Omega.hpp>
+#include<yade/core/Scene.hpp>
+
+YADE_PLUGIN((Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity));
+
+
+
+void Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity::action()// a remplacer par :
+// void Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<InteractionGeometry>& iG, shared_ptr<InteractionPhysics>& iP, Interaction* contact, Scene* scene)
+{
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
+
+ Real dt = Omega::instance().getTimeStep();
+
+ InteractionContainer::iterator ii = scene->interactions->begin(); // a supprimer pr passage au go
+ InteractionContainer::iterator iiEnd = scene->interactions->end(); // a supprimer pr passage au go
+ int nbreInteracTot=0;// a supprimer pr passage au go
+ int nbreInteracMomPlastif=0;// a supprimer pr passage au go
+ for ( ; ii!=iiEnd ; ++ii )// a supprimer pr passage au go
+ {// a supprimer pr passage au go
+ if ((*ii)->interactionGeometry && (*ii)->interactionPhysics)
+ if ((*ii)->isReal())// a supprimer pr passage au go
+ {
+ nbreInteracTot++;
+ const shared_ptr<Interaction>& contact = *ii;// supprimable
+ int id1 = contact->getId1();
+ int id2 = contact->getId2();
+ cout << "contact entre " << id1 << " et " << id2 << " reel ? " << contact->isReal() << endl;
+ if ( !( (*bodies)[id1]->getGroupMask() & (*bodies)[id2]->getGroupMask() & sdecGroupMask) )
+ continue; // skip other groups,
+
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
+ ScGeom* currentContactGeometry = YADE_CAST<ScGeom*>(contact->interactionGeometry.get());
+ NormalInelasticityPhys* currentContactPhysics = YADE_CAST<NormalInelasticityPhys*> (contact->interactionPhysics.get());
+
+ 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;
+ 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
+ currentContactPhysics->normalForce = Vector3r::ZERO;
+ currentContactPhysics->shearForce = Vector3r::ZERO;
+
+ //return;
+ // else
+ // currentContactPhysics->normalForce = -currentContactPhysics->normalAdhesion*currentContactGeometry->normal;
+ }
+ else
+ {
+
+ Vector3r axis;
+ Real angle;
+
+// Here is the code with approximated rotations
+
+ axis = currentContactPhysics->prevNormal.Cross(currentContactGeometry->normal);
+ shearForce -= shearForce.Cross(axis);
+ angle = dt*0.5*currentContactGeometry->normal.Dot(de1->angVel+de2->angVel);
+ axis = angle*currentContactGeometry->normal;
+ shearForce -= shearForce.Cross(axis);
+
+// Here is the code with exact rotationns
+
+ // Quaternionr q;
+ //
+ // axis = currentContactPhysics->prevNormal.cross(currentContactGeometry->normal);
+ // angle = acos(currentContactGeometry->normal.dot(currentContactPhysics->prevNormal));
+ // q.fromAngleAxis(angle,axis);
+ //
+ // currentContactPhysics->shearForce = currentContactPhysics->shearForce*q;
+ //
+ // angle = dt*0.5*currentContactGeometry->normal.dot(de1->angVel+de2->angVel);
+ // axis = currentContactGeometry->normal;
+ // 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)
+ 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_));
+ Vector3r shearVelocity = relativeVelocity-currentContactGeometry->normal.Dot(relativeVelocity)*currentContactGeometry->normal;
+ Vector3r shearDisplacement = shearVelocity*dt;
+
+
+ shearForce -= currentContactPhysics->ks*shearDisplacement;
+/////////////////////// CREEP END /////////////////////////////
+
+ // cerr << "shearForce = " << shearForce << endl;
+ Real maxFs = 0;
+ Real Fs = currentContactPhysics->shearForce.Length();
+ maxFs = std::max((Real) 0,Fn*currentContactPhysics->tangensOfFrictionAngle);
+
+ if ( Fs > maxFs )
+ {
+
+
+ currentContactPhysics->SetBreakingState();
+
+ // maxFs = currentContactPhysics->shearAdhesion;
+ // if (!currentContactPhysics->cohesionDisablesFriction && un>0)
+ // maxFs += Fn * currentContactPhysics->tangensOfFrictionAngle;
+
+ maxFs = max((Real) 0, Fn * currentContactPhysics->tangensOfFrictionAngle);
+
+ //cerr << "currentContactPhysics->tangensOfFrictionAngle = " << currentContactPhysics->tangensOfFrictionAngle << endl;
+ // cerr << "maxFs = " << maxFs << endl;
+
+ maxFs = maxFs / Fs;
+ // cerr << "maxFs = " << maxFs << endl;
+ if (maxFs>1)
+ cerr << "maxFs>1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << endl;
+ shearForce *= maxFs;
+ if (Fn<0) currentContactPhysics->normalForce = Vector3r::ZERO;
+ }
+
+
+
+ //////// PFC3d SlipModel
+
+ Vector3r f = currentContactPhysics->normalForce + shearForce;
+ scene->forces.addForce (id1,-f);
+ scene->forces.addForce (id2,+f);
+ scene->forces.addTorque(id1,-c1x.Cross(f));
+ scene->forces.addTorque(id2, c2x.Cross(f));
+
+
+ // Moment law ///
+ if(momentRotationLaw)
+ {
+ {// updates only orientation of contact (local coordinate system)
+ Vector3r axis = currentContactPhysics->prevNormal.UnitCross(currentContactGeometry->normal);
+ Real angle = unitVectorsAngle(currentContactPhysics->prevNormal,currentContactGeometry->normal);
+ Quaternionr align(axis,angle);
+ currentContactPhysics->currentContactOrientation = align * currentContactPhysics->currentContactOrientation;
+ }
+
+ Quaternionr delta( de1->se3.orientation * currentContactPhysics->initialOrientation1.Conjugate() *
+ currentContactPhysics->initialOrientation2 * de2->se3.orientation.Conjugate());
+
+ Vector3r axis; // axis of rotation - this is the Moment direction UNIT vector.
+ Real angle; // angle represents the power of resistant ELASTIC moment
+ delta.ToAxisAngle(axis,angle);
+ if(angle > Mathr::PI) angle -= Mathr::TWO_PI; // angle is between 0 and 2*pi, but should be between -pi and pi
+
+ //This indentation is a rewrite of original equations (the two commented lines), should work exactly the same.
+//Real elasticMoment = currentContactPhysics->kr * std::abs(angle); // positive value (*)
+
+ Real angle_twist(angle * axis.Dot(currentContactGeometry->normal) );
+ Vector3r axis_twist(angle_twist * currentContactGeometry->normal);
+ Vector3r moment_twist(axis_twist * currentContactPhysics->kr);
+
+ Vector3r axis_bending(angle*axis - axis_twist);
+ Vector3r moment_bending(axis_bending * currentContactPhysics->kr);
+
+//Vector3r moment = axis * elasticMoment * (angle<0.0?-1.0:1.0); // restore sign. (*)
+
+ Vector3r moment = moment_twist + moment_bending;
+
+// Limitation par seuil plastique
+ if (!momentAlwaysElastic)
+ {
+ Real normeMomentMax = currentContactPhysics->forMaxMoment * std::fabs(Fn);
+ Real normeMoment = moment.Length();
+ if(normeMoment>normeMomentMax)
+ {
+ moment *= normeMomentMax/normeMoment;
+ nbreInteracMomPlastif++;
+ }
+ }
+ scene->forces.addTorque(id1,-moment);
+ scene->forces.addTorque(id2,+moment);
+ }
+ // Moment law END ///
+
+ currentContactPhysics->prevNormal = currentContactGeometry->normal;
+// }
+// }
+ }
+// cout << " En tout :"<< nbreInteracTot << "interactions (reelles)" << endl;
+ //cerr << "ncount= " << ncount << endl;//REMOVE
+// cout << momentAlwaysElastic << endl;
+// cout << "Sur " << nbreInteracTot << " interactions (reelles) " << nbreInteracMomPlastif << " se sont vues corriger leur moment" << endl;
+}
+}
+}
=== added file 'pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp 1970-01-01 00:00:00 +0000
+++ pkg/dem/Engine/GlobalEngine/NormalInelasticityLaw.hpp 2010-04-08 09:14:11 +0000
@@ -0,0 +1,37 @@
+/*************************************************************************
+* 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/core/GlobalEngine.hpp> // a remplacer par :
+#include<yade/pkg-common/LawFunctor.hpp>
+#include <set>
+#include <boost/tuple/tuple.hpp>
+
+
+
+class Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity : public GlobalEngine // a remplacer par :
+// class Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity : public LawFunctor
+{
+ public :
+ void action();// a remplacer par :
+// virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, Scene*);
+
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity,
+// LawFunctor,
+ GlobalEngine,
+ "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 (Annual Report 2006) in http://geo.hmg.inpg.fr/frederic/Discrete_Element_Group_FVD.html 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:`CohesiveFrictionalMat` for bodies, with *isCohesive* = 1 (A verifier ce dernier point)\n- :yref:`Ip2_2xCohFrictMat_NormalInelasticityPhys` (=> which involves interactions of :yref:`NormalInelasticityPhys` type)",
+ ((int,sdecGroupMask,1,"?"))
+ ((Real,coeff_dech,1.0,"=kn(unload) / kn(load)"))
+ ((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);
+
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNLEngine.hpp'
--- pkg/dem/Engine/PartialEngine/KinemCNLEngine.hpp 2010-04-03 16:40:33 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNLEngine.hpp 2010-04-08 09:14:11 +0000
@@ -12,7 +12,7 @@
#include<yade/core/PartialEngine.hpp>
#include<yade/core/Body.hpp>
#include<yade/lib-base/Math.hpp>
-#include<yade/pkg-dem/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp>
+#include<yade/pkg-dem/NormalInelasticityLaw.hpp>
=== modified file 'pkg/dem/Engine/PartialEngine/KinemCNSEngine.hpp'
--- pkg/dem/Engine/PartialEngine/KinemCNSEngine.hpp 2010-04-03 16:40:33 +0000
+++ pkg/dem/Engine/PartialEngine/KinemCNSEngine.hpp 2010-04-08 09:14:11 +0000
@@ -12,7 +12,7 @@
#include<yade/core/PartialEngine.hpp>
#include<yade/core/Body.hpp>
#include<yade/lib-base/Math.hpp>
-#include<yade/pkg-dem/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp>
+#include<yade/pkg-dem/NormalInelasticityLaw.hpp>
=== modified file 'pkg/dem/PreProcessor/SimpleShear.cpp'
--- pkg/dem/PreProcessor/SimpleShear.cpp 2010-03-27 22:18:10 +0000
+++ pkg/dem/PreProcessor/SimpleShear.cpp 2010-04-08 09:14:11 +0000
@@ -14,7 +14,7 @@
#include "SimpleShear.hpp"
#include <yade/pkg-dem/CohesiveFrictionalMat.hpp>
-#include <yade/pkg-dem/Law2_ScGeom_NormalInelasticityPhys_NormalInelasticity.hpp>
+#include<yade/pkg-dem/NormalInelasticityLaw.hpp>
#include <yade/pkg-dem/Ip2_2xCohFrictMat_NormalInelasticityPhys.hpp>
#include<yade/pkg-dem/GlobalStiffnessTimeStepper.hpp>
#include <yade/pkg-dem/PositionOrientationRecorder.hpp>
