yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #04525
[Branch ~yade-dev/yade/trunk] Rev 2245: - Fix the plastic dissipation equation in ElasticContactLaw and make plasticDissipation a OpenMPa...
------------------------------------------------------------
revno: 2245
committer: Bruno Chareyre <bchareyre@r1arduina>
branch nick: trunk
timestamp: Thu 2010-05-20 20:39:42 +0200
message:
- Fix the plastic dissipation equation in ElasticContactLaw and make plasticDissipation a OpenMPaccumulator.
- Fix the name scaleDisplacementT(Real multiplier), handle the case force=0 correctly (thanks Vaclav).
- Fix mail adress.
modified:
pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp
pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp
pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp
pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp
pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp
pkg/dem/DataClass/Material/CohFrictMat.cpp
pkg/dem/DataClass/Material/CohFrictMat.hpp
pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp
pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.hpp
pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp
pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.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/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2010-05-18 16:44:30 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.cpp 2010-05-20 18:39:42 +0000
@@ -79,14 +79,14 @@
}
/*! As above : perform slip of the projected contact points. Here, we directly give the multiplier applied on the distance for faster results.
- * The slipped distance is returned.
+ * The plastic displacement (vector) is returned.
*/
-Real Dem3DofGeom_SphereSphere::scaleToDisplacementTMax(Real multiplier){
+Vector3r Dem3DofGeom_SphereSphere::scaleDisplacementT(Real multiplier){
assert(multiplier>=0 && multiplier<=1);
Vector3r p1=contPtInTgPlane1(), p2=contPtInTgPlane2();
Vector3r diff=.5*(multiplier-1)*(p2-p1);
setTgPlanePts(p1-diff,p2+diff);
- return 0;
+ return 2*diff;
}
=== modified file 'pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp'
--- pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2010-05-18 16:29:56 +0000
+++ pkg/dem/DataClass/InteractionGeometry/Dem3DofGeom_SphereSphere.hpp 2010-05-20 18:39:42 +0000
@@ -28,7 +28,7 @@
#endif
}
Real slipToDisplacementTMax(Real displacementTMax);
- Real scaleToDisplacementTMax(Real multiplier);
+ Vector3r scaleDisplacementT(Real multiplier);
/********* end API ***********/
YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(Dem3DofGeom_SphereSphere,Dem3DofGeom,"Class representing 2 spheres in contact which computes 3 degrees of freedom (normal and shear deformation).",
=== modified file 'pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2010-05-18 16:29:56 +0000
+++ pkg/dem/DataClass/InteractionGeometry/DemXDofGeom.hpp 2010-05-20 18:39:42 +0000
@@ -30,7 +30,7 @@
virtual Real displacementN();
virtual Vector3r displacementT(){throw;}
virtual Real slipToDisplacementTMax(Real displacementTMax){throw;}; // plasticity
- virtual Real scaleToDisplacementTMax(Real multiplier){throw;}; // plasticity (variant using multiplier dispMax/disp)
+ virtual Vector3r scaleDisplacementT(Real multiplier){throw;}; // plasticity (variant using multiplier dispMax/disp)
// reference radii, for contact stiffness computation; set to negative for nonsense values
// end API
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-04-25 15:46:26 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-05-20 18:39:42 +0000
@@ -1,6 +1,7 @@
// © 2004 Olivier Galizzi <olivier.galizzi@xxxxxxx>
// © 2004 Janek Kozicki <cosurgi@xxxxxxxxxx>
// © 2008 Václav Šmilauer <eudoxos@xxxxxxxx>
+// © 2008 Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
#include "ScGeom.hpp"
#include<yade/core/Omega.hpp>
@@ -8,7 +9,6 @@
// At least one virtual method must be in the .cpp file (!!!)
ScGeom::~ScGeom(){};
-#ifdef SCG_SHEAR
Vector3r ScGeom::updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
Vector3r axis;
@@ -63,20 +63,16 @@
shear+=shearIncrement;
return shearIncrement;
}
-#endif
Vector3r ScGeom::updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
-
Vector3r axis;
Real angle;
-
// approximated rotations
axis = prevNormal.cross(normal);
shearForce -= shearForce.cross(axis);
angle = dt*0.5*normal.dot(rbp1->angVel + rbp2->angVel);
axis = angle*normal;
shearForce -= shearForce.cross(axis);
-
// exact rotations
#if 0
Quaternionr q;
@@ -92,7 +88,6 @@
Vector3r& x = contactPoint;
Vector3r c1x, c2x;
-
if(avoidGranularRatcheting){
/* The following definition of c1x and c2x is to avoid "granular ratcheting"
* (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
@@ -108,7 +103,6 @@
c1x = (x - rbp1->pos);
c2x = (x - rbp2->pos);
}
-
Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
Vector3r shearVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
Vector3r shearDisplacement = shearVelocity*dt;
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-05-20 18:39:42 +0000
@@ -20,21 +20,17 @@
Real penetrationDepth;
Real &radius1, &radius2;
- #ifdef SCG_SHEAR
- //! update shear on this contact given dynamic parameters of bodies. Should be called from constitutive law, exactly once per iteration. Returns shear increment in this update, which is already added to shear.
- Vector3r updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
- #endif
+ //! update shear on this contact given dynamic parameters of bodies. Should be called from constitutive law, exactly once per iteration. Returns shear increment in this update, which is already added to shear.
+ Vector3r updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
virtual ~ScGeom();
-
+ //! update the shear force on this contact given dynamic parameters of bodies (rotate the force and add the force increment). Returns shear displacement increment, don't update shear.
Vector3r updateShearForce(Vector3r& shearForce, Real ks, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
- YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`spheres<Sphere>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`spheres<Sphere>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
((Vector3r,contactPoint,Vector3r::Zero(),"Reference point of the contact. |ycomp|"))
- #ifdef SCG_SHEAR
- ((Vector3r,shear,Vector3r::Zero(),"Total value of the current shear. Update the value using ScGeom::updateShear. |ycomp|"))
- ((Vector3r,prevNormal,Vector3r::Zero(),"Normal of the contact in the previous step. |ycomp|"))
- #endif
+ ((Vector3r,shear,Vector3r::Zero(),"Total value of the current shear. Update the value using ScGeom::updateShear. |ycomp|"))
+ ((Vector3r,prevNormal,Vector3r::Zero(),"Normal of the contact in the previous step. |ycomp|"))
,
/* extra initializers */ ((radius1,GenericSpheresContact::refR1)) ((radius2,GenericSpheresContact::refR2)),
/* ctor */ createIndex();,
=== modified file 'pkg/dem/DataClass/Material/CohFrictMat.cpp'
--- pkg/dem/DataClass/Material/CohFrictMat.cpp 2010-04-19 13:08:27 +0000
+++ pkg/dem/DataClass/Material/CohFrictMat.cpp 2010-05-20 18:39:42 +0000
@@ -1,5 +1,5 @@
/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxx> *
+* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
* *
* This program is free software; it is licensed under the terms of the *
=== modified file 'pkg/dem/DataClass/Material/CohFrictMat.hpp'
--- pkg/dem/DataClass/Material/CohFrictMat.hpp 2010-04-19 13:08:27 +0000
+++ pkg/dem/DataClass/Material/CohFrictMat.hpp 2010-05-20 18:39:42 +0000
@@ -1,5 +1,5 @@
/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxx> *
+* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
* *
* This program is free software; it is licensed under the terms of the *
=== modified file 'pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp 2010-05-03 12:17:44 +0000
+++ pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.cpp 2010-05-20 18:39:42 +0000
@@ -1,5 +1,5 @@
/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxx> *
+* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
* *
* This program is free software; it is licensed under the terms of the *
=== modified file 'pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.hpp'
--- pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.hpp 2010-04-19 13:08:27 +0000
+++ pkg/dem/Engine/Functor/Ip2_2xCohFrictMat_CohFrictPhys.hpp 2010-05-20 18:39:42 +0000
@@ -1,5 +1,5 @@
/*************************************************************************
-* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxx> *
+* Copyright (C) 2007 by Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx> *
* Copyright (C) 2008 by Janek Kozicki <cosurgi@xxxxxxxxxx> *
* *
* This program is free software; it is licensed under the terms of the *
=== modified file 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-04-26 09:02:51 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.cpp 2010-05-20 18:39:42 +0000
@@ -31,15 +31,15 @@
Real Va = mat1->poisson;
Real Vb = mat2->poisson;
- Real Da,Db; Vector3r normal;
+ Real Ra,Rb; Vector3r normal;
assert(dynamic_cast<GenericSpheresContact*>(interaction->interactionGeometry.get()));//only in debug mode
GenericSpheresContact* sphCont=YADE_CAST<GenericSpheresContact*>(interaction->interactionGeometry.get());
- {Da=sphCont->refR1>0?sphCont->refR1:sphCont->refR2; Db=sphCont->refR2>0?sphCont->refR2:sphCont->refR1; normal=sphCont->normal;}
+ {Ra=sphCont->refR1>0?sphCont->refR1:sphCont->refR2; Rb=sphCont->refR2>0?sphCont->refR2:sphCont->refR1; normal=sphCont->normal;}
- //harmonic average of the two stiffnesses when (Di.Ei/2) is the stiffness of a contact point on sphere "i"
- Real Kn = 2*Ea*Da*Eb*Db/(Ea*Da+Eb*Db);
+ //harmonic average of the two stiffnesses when (Ri.Ei/2) is the stiffness of a contact point on sphere "i"
+ Real Kn = 2*Ea*Ra*Eb*Rb/(Ea*Ra+Eb*Rb);
//same for shear stiffness
- Real Ks = 2*Ea*Da*Va*Eb*Db*Vb/(Ea*Da*Va+Eb*Db*Va);
+ Real Ks = 2*Ea*Ra*Va*Eb*Rb*Vb/(Ea*Ra*Va+Eb*Rb*Va);
contactPhysics->frictionAngle = std::min(mat1->frictionAngle,mat2->frictionAngle);
contactPhysics->tangensOfFrictionAngle = std::tan(contactPhysics->frictionAngle);
=== modified file 'pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp'
--- pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp 2010-05-02 16:02:29 +0000
+++ pkg/dem/Engine/Functor/Ip2_FrictMat_FrictMat_FrictPhys.hpp 2010-05-20 18:39:42 +0000
@@ -16,7 +16,7 @@
const shared_ptr<Material>& b2,
const shared_ptr<Interaction>& interaction);
FUNCTOR2D(FrictMat,FrictMat);
- YADE_CLASS_BASE_DOC(Ip2_FrictMat_FrictMat_FrictPhys,InteractionPhysicsFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under symetric point loads is defined here as 1/(E.D), with E the stiffness of the sphere and D its diameter, and corresponds to a compliance 1/(2.E.D) from each contact point. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. 1/(2.E.D1)+1/(2.E.D1) in the general case, or 1/(E.D) in the special case of equal sizes. Note that summing compliances corresponds to an harmonic average of stiffnesss, which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor. The friction angle of the contact is defined as the minimum angle of the two materials in contact.\n\nOnly interactions with :yref:`ScGeom` or :yref:`Dem3DofGeom` geometry are meaningfully accepted; run-time typecheck can make this functor unnecessarily slow in general. Such design is problematic in itself, though -- from http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg02603.html:\n\n\tYou have to suppose some exact type of InteractionGeometry in the Ip2 functor, but you don't know anything about it (Ip2 only guarantees you get certain InteractionPhysics types, via the dispatch mechanism).\n\n\tThat means, unless you use Ig2 functor producing the desired type, the code will break (crash or whatever). The right behavior would be either to accept any type (what we have now, at least in principle), or really enforce InteractionGeometry type of the interation passed to that particular Ip2 functor.\n\nEtc.");
+ YADE_CLASS_BASE_DOC(Ip2_FrictMat_FrictMat_FrictPhys,InteractionPhysicsFunctor,"Create a :yref:`FrictPhys` from two :yref:`FrictMats<FrictMat>`. The compliance of one sphere under symetric point loads is defined here as 1/(E.r), with E the stiffness of the sphere and r its radius, and corresponds to a compliance 1/(2.E.r)=1/(E.D) from each contact point. The compliance of the contact itself will be the sum of compliances from each sphere, i.e. 1/(E.D1)+1/(E.D2) in the general case, or 1/(E.r) in the special case of equal sizes. Note that summing compliances corresponds to an harmonic average of stiffnesss, which is how kn is actually computed in the :yref:`Ip2_FrictMat_FrictMat_FrictPhys` functor.\n\n The shear stiffness ks of one sphere is defined via the material parameter :yref:`FrictPhys::poisson`, as ks=poisson*kn, and the resulting shear stiffness of the interaction will be also an harmonic average.\n\n The friction angle of the contact is defined as the minimum angle of the two materials in contact.\n\nOnly interactions with :yref:`ScGeom` or :yref:`Dem3DofGeom` geometry are meaningfully accepted; run-time typecheck can make this functor unnecessarily slow in general. Such design is problematic in itself, though -- from http://www.mail-archive.com/yade-dev@xxxxxxxxxxxxxxxxxxx/msg02603.html:\n\n\tYou have to suppose some exact type of InteractionGeometry in the Ip2 functor, but you don't know anything about it (Ip2 only guarantees you get certain InteractionPhysics types, via the dispatch mechanism).\n\n\tThat means, unless you use Ig2 functor producing the desired type, the code will break (crash or whatever). The right behavior would be either to accept any type (what we have now, at least in principle), or really enforce InteractionGeometry type of the interation passed to that particular Ip2 functor.\n\nEtc.");
};
REGISTER_SERIALIZABLE(Ip2_FrictMat_FrictMat_FrictPhys);
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-05-18 16:29:56 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-05-20 18:39:42 +0000
@@ -1,6 +1,6 @@
/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi *
-* olivier.galizzi@xxxxxxx *
+* Copyright (C) 2004 by Olivier Galizzi olivier.galizzi@xxxxxxx *
+* Copyright (C) 2009 by Bruno Chareyre bruno.chareyre@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. *
@@ -16,6 +16,9 @@
YADE_PLUGIN((Law2_ScGeom_FrictPhys_Basic)(Law2_Dem3DofGeom_FrictPhys_Basic)(ElasticContactLaw)(Law2_Dem6DofGeom_FrictPhys_Beam));
+Real Law2_ScGeom_FrictPhys_Basic::Real0=0;
+Real Law2_ScGeom_FrictPhys_Basic::getPlasticDissipation() {return (Real) plasticDissipation;}
+void Law2_ScGeom_FrictPhys_Basic::initPlasticDissipation(Real initVal) {plasticDissipation.reset(); plasticDissipation+=initVal;}
void ElasticContactLaw::action()
{
@@ -64,7 +67,8 @@
TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
- if (!traceEnergy){//Update force but don't compute energy terms
+ if (!traceEnergy){
+ //Update force but don't compute energy terms (see below))
if(useShear){
currentContactGeometry->updateShear(de1,de2,dt);
shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
@@ -77,26 +81,21 @@
Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
shearForce *= ratio;
if(useShear) currentContactGeometry->shear*=ratio;}
- } else {//almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false
+ } else {
+ //almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false
+ if(useShear) throw ("energy tracing not defined with useShear==true");
Vector3r prevForce=shearForce;//store prev force for definition of plastic slip
- if(useShear) throw ("energy tracing not defined with useShear==true");
- /*{
- currentContactGeometry->updateShear(de1,de2,dt);
- shearForce=currentContactPhysics->ks*currentContactGeometry->shear;//FIXME : energy terms if useShear?
- } else {*/
Vector3r shearDisp = currentContactGeometry->updateShearForce(shearForce,currentContactPhysics->ks,currentContactPhysics->prevNormal,de1,de2,dt);
- //}
- // PFC3d SlipModel, is using friction angle. CoulombCriterion
Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
if( shearForce.squaredNorm() > maxFs ){
Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
+ Vector3r trialForce=shearForce;//store prev force for definition of plastic slip
//define the plastic work input and increment the total plastic energy dissipated
+ shearForce *= ratio;
plasticDissipation +=
- (shearDisp+(1/currentContactPhysics->ks)*(shearForce-prevForce))//plastic disp.
+ ((1/currentContactPhysics->ks)*(trialForce-shearForce))//plastic disp.
.dot(shearForce);//active force
- shearForce *= ratio;
- //if(useShear) currentContactGeometry->shear*=ratio;
}
}
applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
@@ -112,12 +111,10 @@
phys->normalForce=phys->kn*displN*geom->normal;
Real maxFsSq=phys->normalForce.squaredNorm()*pow(phys->tangensOfFrictionAngle,2);
Vector3r trialFs=phys->ks*geom->displacementT();
- Real multiplier=maxFsSq/trialFs.squaredNorm();
-// if(trialFs.squaredNorm()>maxFsSq){
-// geom->slipToDisplacementTMax(sqrt(maxFsSq)/phys->ks); trialFs*=sqrt(maxFsSq/(trialFs.squaredNorm()));}
- if(multiplier<1){
- multiplier = sqrt(multiplier);
- geom->scaleToDisplacementTMax(multiplier); trialFs*=multiplier;}
+ Real trialFsSq = trialFs.squaredNorm();
+ if(trialFsSq>maxFsSq){
+ Real multiplier=sqrt(maxFsSq/trialFsSq);
+ geom->scaleDisplacementT(multiplier); trialFs*=multiplier;}
phys->shearForce=trialFs;
applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,scene);
}
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-05-18 10:28:02 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-05-20 18:39:42 +0000
@@ -1,6 +1,6 @@
/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi *
-* olivier.galizzi@xxxxxxx *
+* Copyright (C) 2004 by Olivier Galizzi olivier.galizzi@xxxxxxx *
+* Copyright (C) 2009 by Bruno Chareyre bruno.chareyre@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. *
@@ -9,26 +9,28 @@
#pragma once
#include<yade/core/GlobalEngine.hpp>
-
-// only to see whether SCG_SHEAR is defined, may be removed in the future
-#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/LawFunctor.hpp>
#include <set>
#include <boost/tuple/tuple.hpp>
-
+#include<yade/lib-base/openmp-accu.hpp>
class Law2_ScGeom_FrictPhys_Basic: public LawFunctor{
public:
+ static Real Real0;
+ OpenMPAccumulator<Real,&Law2_ScGeom_FrictPhys_Basic::Real0> plasticDissipation;
virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
Real elasticEnergy ();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom`.",
+ Real getPlasticDissipation();
+ void initPlasticDissipation(Real initVal=0);
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the latest).",
((bool,neverErase,false,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
((bool,useShear,false,"Use ScGeom::updateShear rather than ScGeom::updateShearForce for shear force computation."))
((bool,traceEnergy,false,"Define the total energy dissipated in plastic slips at all contacts."))
- ((Real,plasticDissipation,0,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_Basic::traceEnergy` is true. |yupdate|"))
,,
.def("elasticEnergy",&Law2_ScGeom_FrictPhys_Basic::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
+ .def("plasticDissipation",&Law2_ScGeom_FrictPhys_Basic::getPlasticDissipation,"Total energy dissipated in plastic slips at all FrictPhys contacts. Computed only if :yref:`Law2_ScGeom_FrictPhys_Basic::traceEnergy` is true.")
+ .def("initPlasticDissipation",&Law2_ScGeom_FrictPhys_Basic::initPlasticDissipation,"Initialize cummulated plastic dissipation to a value (0 by default).")
);
FUNCTOR2D(ScGeom,FrictPhys);
DECLARE_LOGGER;
