yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #11085
[Branch ~yade-pkg/yade/git-trunk] Rev 4091: Law2 return bool - fix https://bugs.launchpad.net/yade/+bug/1324190
------------------------------------------------------------
revno: 4091
committer: Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
timestamp: Fri 2014-07-18 20:18:50 +0200
message:
Law2 return bool - fix https://bugs.launchpad.net/yade/+bug/1324190
modified:
pkg/common/Cylinder.cpp
pkg/common/Cylinder.hpp
pkg/common/Dispatching.hpp
pkg/common/Grid.cpp
pkg/common/Grid.hpp
pkg/common/InteractionLoop.cpp
pkg/dem/BubbleMat.cpp
pkg/dem/BubbleMat.hpp
pkg/dem/CohesiveFrictionalContactLaw.cpp
pkg/dem/CohesiveFrictionalContactLaw.hpp
pkg/dem/ConcretePM.cpp
pkg/dem/ConcretePM.hpp
pkg/dem/ElasticContactLaw.cpp
pkg/dem/ElasticContactLaw.hpp
pkg/dem/FrictViscoPM.cpp
pkg/dem/FrictViscoPM.hpp
pkg/dem/HertzMindlin.cpp
pkg/dem/HertzMindlin.hpp
pkg/dem/InelastCohFrictPM.cpp
pkg/dem/InelastCohFrictPM.hpp
pkg/dem/JointedCohesiveFrictionalPM.cpp
pkg/dem/JointedCohesiveFrictionalPM.hpp
pkg/dem/L3Geom.cpp
pkg/dem/L3Geom.hpp
pkg/dem/LudingPM.cpp
pkg/dem/LudingPM.hpp
pkg/dem/NormalInelasticPM.cpp
pkg/dem/NormalInelasticPM.hpp
pkg/dem/Polyhedra.cpp
pkg/dem/Polyhedra.hpp
pkg/dem/Tetra.cpp
pkg/dem/Tetra.hpp
pkg/dem/ViscoelasticCapillarPM.cpp
pkg/dem/ViscoelasticCapillarPM.hpp
pkg/dem/ViscoelasticPM.cpp
pkg/dem/ViscoelasticPM.hpp
pkg/dem/WirePM.cpp
pkg/dem/WirePM.hpp
--
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
=== modified file 'pkg/common/Cylinder.cpp'
--- pkg/common/Cylinder.cpp 2014-07-03 07:16:58 +0000
+++ pkg/common/Cylinder.cpp 2014-07-18 18:18:50 +0000
@@ -658,7 +658,7 @@
}
}
-void Law2_CylScGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_CylScGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
CylScGeom* geom= static_cast<CylScGeom*>(ig.get());
@@ -667,13 +667,11 @@
if (neverErase) {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();}
- else scene->interactions->requestErase(contact);
- return;}
+ else return false;}
if (geom->isDuplicate) {
if (id2!=geom->trueInt) {
//cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
- if (geom->isDuplicate==2) {/*cerr<<"erase duplicate "<<id1<<" "<<id2<<endl;*/scene->interactions->requestErase(contact);}
- return;}
+ if (geom->isDuplicate==2) return false;}
}
Real& un=geom->penetrationDepth;
phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
@@ -721,10 +719,11 @@
scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
}
+ return true;
}
-void Law2_CylScGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact) {
+bool Law2_CylScGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact) {
int id1 = contact->getId1(), id2 = contact->getId2();
@@ -738,19 +737,17 @@
if (geom->isDuplicate) {
if (id2!=geom->trueInt) {
//cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
- if (geom->isDuplicate==2) {/*cerr<<"erase duplicate coh "<<id1<<" "<<id2<<endl;*/scene->interactions->requestErase(contact);}
- return;}
+ if (geom->isDuplicate==2) return false;}
}
if (currentContactPhysics->fragile && (-Fn)> currentContactPhysics->normalAdhesion) {
// BREAK due to tension
- scene->interactions->requestErase(contact);
+ return false;
} else {
if ((-Fn)> currentContactPhysics->normalAdhesion) {//normal plasticity
Fn=-currentContactPhysics->normalAdhesion;
currentContactPhysics->unp = un+currentContactPhysics->normalAdhesion/currentContactPhysics->kn;
- if (currentContactPhysics->unpMax && currentContactPhysics->unp<currentContactPhysics->unpMax)
- scene->interactions->requestErase(contact);
+ if (currentContactPhysics->unpMax && currentContactPhysics->unp<currentContactPhysics->unpMax) return false;
}
currentContactPhysics->normalForce = Fn*geom->normal;
Vector3r& shearForce = geom->rotate(currentContactPhysics->shearForce);
@@ -794,11 +791,12 @@
}
//applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
}
+ return true;
}
-void Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
ChCylGeom6D* geom= YADE_CAST<ChCylGeom6D*>(ig.get());
@@ -816,24 +814,15 @@
if (contact->isFresh(scene)) shearForce = Vector3r::Zero(); //contact nouveau => force tengentielle = 0,0,0
Real un = geom->penetrationDepth; //un : interpenetration
Real Fn = currentContactPhysics->kn*(un-currentContactPhysics->unp); //Fn : force normale
- /*if (geom->isDuplicate) {
- if (id2!=geom->trueInt) {
- //cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
- if (geom->isDuplicate==2) {cerr<<"erase duplicate coh "<<id1<<" "<<id2<<endl;scene->interactions->requestErase(contact);}
- return;}
- }
-*/
- if (currentContactPhysics->fragile && (-Fn)> currentContactPhysics->normalAdhesion) {
- // BREAK due to tension
- scene->interactions->requestErase(contact);
- } else {
+ if (currentContactPhysics->fragile && (-Fn)> currentContactPhysics->normalAdhesion) return false; // BREAK due to tension
+ else {
if ((-Fn)> currentContactPhysics->normalAdhesion) {//normal plasticity
Fn=-currentContactPhysics->normalAdhesion;
currentContactPhysics->unp = un+currentContactPhysics->normalAdhesion/currentContactPhysics->kn;
if (currentContactPhysics->unpMax && currentContactPhysics->unp<currentContactPhysics->unpMax)
- scene->interactions->requestErase(contact);
- }
+ return false;
+ }
currentContactPhysics->normalForce = Fn*geom->normal;
@@ -885,6 +874,7 @@
scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
}
}
+ return true;
}
=== modified file 'pkg/common/Cylinder.hpp'
--- pkg/common/Cylinder.hpp 2014-07-14 20:08:33 +0000
+++ pkg/common/Cylinder.hpp 2014-07-18 18:18:50 +0000
@@ -218,7 +218,7 @@
class Law2_CylScGeom_FrictPhys_CundallStrack: public LawFunctor{
public:
//OpenMPAccumulator<Real> plasticDissipation;
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
//Real elasticEnergy ();
//Real getPlasticDissipation();
//void initPlasticDissipation(Real initVal=0);
@@ -239,7 +239,7 @@
class Law2_CylScGeom6D_CohFrictPhys_CohesionMoment: public LawFunctor {
public:
//OpenMPAccumulator<Real> plasticDissipation;
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
//Real elasticEnergy ();
//Real getPlasticDissipation();
//void initPlasticDissipation(Real initVal=0);
@@ -265,7 +265,7 @@
class Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment: public LawFunctor {
public:
//OpenMPAccumulator<Real> plasticDissipation;
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ChCylGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear compression, and Mohr-Coulomb plasticity surface without cohesion.\nThis law implements the classical linear elastic-plastic law from [CundallStrack1979]_ (see also [Pfc3dManual30]_). The normal force is (with the convention of positive tensile forces) $F_n=\\min(k_n u_n, 0)$. The shear force is $F_s=k_s u_s$, the plasticity condition defines the maximum value of the shear force : $F_s^{\\max}=F_n\\tan(\\phi)$, with $\\phi$ the friction angle.\n\n.. note::\n This law is well tested in the context of triaxial simulation, and has been used for a number of published results (see e.g. [Scholtes2009b]_ and other papers from the same authors). It is generalised by :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`, which adds cohesion and moments at contact.",
((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,traceEnergy,false,Attr::hidden,"Define the total energy dissipated in plastic slips at all contacts."))
=== modified file 'pkg/common/Dispatching.hpp'
--- pkg/common/Dispatching.hpp 2014-02-22 04:56:33 +0000
+++ pkg/common/Dispatching.hpp 2014-07-18 18:18:50 +0000
@@ -52,7 +52,7 @@
class LawFunctor: public Functor2D<
/*dispatch types*/ IGeom,IPhys,
- /*return type*/ void,
+ /*return type*/ bool,
/*argument types*/ TYPELIST_3(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*)
>{
public: virtual ~LawFunctor();
=== modified file 'pkg/common/Grid.cpp'
--- pkg/common/Grid.cpp 2014-07-03 11:31:29 +0000
+++ pkg/common/Grid.cpp 2014-07-18 18:18:50 +0000
@@ -397,7 +397,7 @@
//!################## Laws #####################
//! O/
-void Law2_ScGridCoGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGridCoGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
ScGridCoGeom* geom= static_cast<ScGridCoGeom*>(ig.get());
FrictPhys* phys = static_cast<FrictPhys*>(ip.get());
@@ -405,16 +405,11 @@
if (neverErase) {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();}
- else scene->interactions->requestErase(contact);
- return;}
+ else return false;}
if (geom->isDuplicate) {
if (id2!=geom->trueInt) {
//cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
- if (geom->isDuplicate==2) {
- //cerr<<"erase duplicate "<<id1<<" "<<id2<<endl;
- scene->interactions->requestErase(contact);
- }
- return;
+ if (geom->isDuplicate==2) return false;
}
}
Real& un=geom->penetrationDepth;
@@ -455,7 +450,7 @@
YADE_PLUGIN((Law2_ScGridCoGeom_FrictPhys_CundallStrack));
-void Law2_ScGridCoGeom_CohFrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGridCoGeom_CohFrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
const int &id1 = contact->getId1();
const int &id2 = contact->getId2();
ScGridCoGeom* geom = YADE_CAST<ScGridCoGeom*> (ig.get());
@@ -464,11 +459,7 @@
if (geom->isDuplicate) {
if (id2!=geom->trueInt) {
//cerr<<"skip duplicate "<<id1<<" "<<id2<<endl;
- if (geom->isDuplicate==2) {
- //cerr<<"erase duplicate "<<id1<<" "<<id2<<endl;
- scene->interactions->requestErase(contact);
- }
- return;
+ if (geom->isDuplicate==2) return false;
}
}
@@ -479,13 +470,13 @@
if (phys->fragile && (-Fn)> phys->normalAdhesion) {
// BREAK due to tension
- scene->interactions->requestErase(contact); return;
+ return false;
} else {
if ((-Fn)> phys->normalAdhesion) {//normal plasticity
Fn=-phys->normalAdhesion;
phys->unp = un+phys->normalAdhesion/phys->kn;
if (phys->unpMax && phys->unp<phys->unpMax)
- scene->interactions->requestErase(contact); return;
+ return false;
}
phys->normalForce = Fn*geom->normal;
Vector3r& shearForce = geom->rotate(phys->shearForce);
@@ -524,7 +515,7 @@
}
YADE_PLUGIN((Law2_ScGridCoGeom_CohFrictPhys_CundallStrack));
-void Law2_GridCoGridCoGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_GridCoGridCoGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
id_t id11 = (static_cast<GridConnection*>((&Body::byId(id1)->shape)->get()))->node1->getId();
id_t id12 = (static_cast<GridConnection*>((&Body::byId(id1)->shape)->get()))->node2->getId();
@@ -536,8 +527,7 @@
if (neverErase) {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();}
- else scene->interactions->requestErase(contact);
- return;}
+ else return false;}
Real& un=geom->penetrationDepth;
phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
@@ -579,6 +569,7 @@
scene->forces.addTorque(id12,geom->relPos1*torque1);
scene->forces.addTorque(id21,(1-geom->relPos2)*torque2);
scene->forces.addTorque(id22,geom->relPos2*torque2);
+ return true;
}
YADE_PLUGIN((Law2_GridCoGridCoGeom_FrictPhys_CundallStrack));
//!################## Bounds #####################
=== modified file 'pkg/common/Grid.hpp'
--- pkg/common/Grid.hpp 2014-07-14 20:08:33 +0000
+++ pkg/common/Grid.hpp 2014-07-18 18:18:50 +0000
@@ -168,7 +168,7 @@
//! O/
class Law2_ScGridCoGeom_FrictPhys_CundallStrack: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGridCoGeom_FrictPhys_CundallStrack,LawFunctor,"Law between a frictional :yref:`GridConnection` and a frictional :yref:`Sphere`. Almost the same than :yref:`Law2_ScGeom_FrictPhys_CundallStrack`, but the force is divided and applied on the two :yref:`GridNodes<GridNode>` only.",
((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`)"))
((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
@@ -181,7 +181,7 @@
class Law2_ScGridCoGeom_CohFrictPhys_CundallStrack: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGridCoGeom_CohFrictPhys_CundallStrack,LawFunctor,"Law between a cohesive frictional :yref:`GridConnection` and a cohesive frictional :yref:`Sphere`. Almost the same than :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment`, but THE ROTATIONAL MOMENTS ARE NOT COMPUTED.",
((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`)"))
((int,plastDissipIx,-1,(Attr::hidden|Attr::noSave),"Index for plastic dissipation (with O.trackEnergy)"))
@@ -194,7 +194,7 @@
//! -/-
class Law2_GridCoGridCoGeom_FrictPhys_CundallStrack: public Law2_ScGeom_FrictPhys_CundallStrack{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS(Law2_GridCoGridCoGeom_FrictPhys_CundallStrack,Law2_ScGeom_FrictPhys_CundallStrack,"Frictional elastic contact law between two :yref:`gridConnection` . See :yref:`Law2_ScGeom_FrictPhys_CundallStrack` for more details.",
/*ATTRS*/
);
=== modified file 'pkg/common/InteractionLoop.cpp'
--- pkg/common/InteractionLoop.cpp 2014-07-02 16:11:24 +0000
+++ pkg/common/InteractionLoop.cpp 2014-07-18 18:18:50 +0000
@@ -137,7 +137,8 @@
assert(!swap); // reverse call would make no sense, as the arguments are of different types
}
assert(I->functorCache.constLaw);
- I->functorCache.constLaw->go(I->geom,I->phys,I.get());
+ //If the functor return false, the interaction is reset
+ if (!I->functorCache.constLaw->go(I->geom,I->phys,I.get())) scene->interactions->requestErase(I);
// process callbacks for this interaction
if(!I->isReal()) continue; // it is possible that Law2_ functor called requestErase, hence this check
=== modified file 'pkg/dem/BubbleMat.cpp'
--- pkg/dem/BubbleMat.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/BubbleMat.cpp 2014-07-18 18:18:50 +0000
@@ -51,7 +51,7 @@
/********************** Law2_ScGeom_BubblePhys_Bubble ****************************/
CREATE_LOGGER(Law2_ScGeom_BubblePhys_Bubble);
-void Law2_ScGeom_BubblePhys_Bubble::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
+bool Law2_ScGeom_BubblePhys_Bubble::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
ScGeom* geom=static_cast<ScGeom*>(_geom.get());
BubblePhys* phys=static_cast<BubblePhys*>(_phys.get());
=== modified file 'pkg/dem/BubbleMat.hpp'
--- pkg/dem/BubbleMat.hpp 2013-07-05 15:54:47 +0000
+++ pkg/dem/BubbleMat.hpp 2014-07-18 18:18:50 +0000
@@ -64,7 +64,7 @@
/********************** Law2_ScGeom_BubblePhys_Bubble ****************************/
class Law2_ScGeom_BubblePhys_Bubble : public LawFunctor{
public:
- void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* interaction);
+ bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* interaction);
FUNCTOR2D(GenericSpheresContact,BubblePhys);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_BubblePhys_Bubble,LawFunctor,"Constitutive law for Bubble model.",
,
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/CohesiveFrictionalContactLaw.cpp 2014-07-14 20:08:33 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.cpp 2014-07-18 18:18:50 +0000
@@ -99,7 +99,7 @@
functor->go(I->geom, I->phys, I.get());}
}
-void Law2_ScGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
+bool Law2_ScGeom6D_CohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
{
const Real& dt = scene->dt;
const int &id1 = contact->getId1();
@@ -114,13 +114,13 @@
if (phys->fragile && (-Fn)> phys->normalAdhesion) {
// BREAK due to tension
- scene->interactions->requestErase(contact); return;
+ return false;
} else {
if ((-Fn)> phys->normalAdhesion) {//normal plasticity
Fn=-phys->normalAdhesion;
phys->unp = un+phys->normalAdhesion/phys->kn;
if (phys->unpMax && phys->unp<phys->unpMax)
- scene->interactions->requestErase(contact); return;
+ return false;
}
phys->normalForce = Fn*geom->normal;
State* de1 = Body::byId(id1,scene)->state.get();
@@ -223,6 +223,7 @@
}
/// Moment law END ///
}
+ return true;
}
=== modified file 'pkg/dem/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-07-14 20:08:33 +0000
+++ pkg/dem/CohesiveFrictionalContactLaw.hpp 2014-07-18 18:18:50 +0000
@@ -86,7 +86,7 @@
Real totalElastEnergy();
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_CohFrictPhys_CohesionMoment,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. This law adds adhesion and moments to :yref:`Law2_ScGeom_FrictPhys_CundallStrack`.\n\nThe normal force is (with the convention of positive tensile forces) $F_n=min(k_n*u_n, a_n)$, with $a_n$ the normal adhesion. The shear force is $F_s=k_s*u_s$, the plasticity condition defines the maximum value of the shear force, by default $F_s^{max}=F_n*tan(\\phi)+a_s$, with $\\phi$ the friction angle and $a_s$ the shear adhesion. If :yref:`CohFrictPhys::cohesionDisableFriction` is True, friction is ignored as long as adhesion is active, and the maximum shear force is only $F_s^{max}=a_s$.\n\nIf the maximum tensile or maximum shear force is reached and :yref:`CohFrictPhys::fragile` =True (default), the cohesive link is broken, and $a_n, a_s$ are set back to zero. If a tensile force is present, the contact is lost, else the shear strength is $F_s^{max}=F_n*tan(\\phi)$. If :yref:`CohFrictPhys::fragile` =False (in course of implementation), the behaviour is perfectly plastic, and the shear strength is kept constant.\n\nIf :yref:`Law2_ScGeom6D_CohFrictPhys_CohesionMoment::momentRotationLaw` =True, bending and twisting moments are computed using a linear law with moduli respectively $k_t$ and $k_r$ (the two values are the same currently), so that the moments are : $M_b=k_b*\\Theta_b$ and $M_t=k_t*\\Theta_t$, with $\\Theta_{b,t}$ the relative rotations between interacting bodies. The maximum value of moments can be defined and takes the form of rolling friction. Cohesive -type moment may also be included in the future.\n\nCreep at contact is implemented in this law, as defined in [Hassan2010]_. If activated, there is a viscous behaviour of the shear and twisting components, and the evolution of the elastic parts of shear displacement and relative twist is given by $du_{s,e}/dt=-F_s/\\nu_s$ and $d\\Theta_{t,e}/dt=-M_t/\\nu_t$.",
((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,always_use_moment_law,false,,"If true, use bending/twisting moments at all contacts. If false, compute moments only for cohesive contacts."))
=== modified file 'pkg/dem/ConcretePM.cpp'
--- pkg/dem/ConcretePM.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/ConcretePM.cpp 2014-07-18 18:18:50 +0000
@@ -285,7 +285,7 @@
#define NNAN(a) YADE_VERIFY(!isnan(a));
#define NNANV(v) YADE_VERIFY(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
-void Law2_ScGeom_CpmPhys_Cpm::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
+bool Law2_ScGeom_CpmPhys_Cpm::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
TIMING_DELTAS_START();
ScGeom* geom=static_cast<ScGeom*>(_geom.get());
CpmPhys* phys=static_cast<CpmPhys*>(_phys.get());
@@ -413,8 +413,7 @@
{ boost::mutex::scoped_lock lock(st1->updateMutex); st1->numBrokenCohesive += 1; /* st1->epsPlBroken += epsPlSum; */ }
{ boost::mutex::scoped_lock lock(st2->updateMutex); st2->numBrokenCohesive += 1; /* st2->epsPlBroken += epsPlSum; */ }
/* } */
- scene->interactions->requestErase(I);
- return;
+ return false;
}
Fn = sigmaN*crossSection; phys->normalForce = -Fn*geom->normal;
@@ -438,6 +437,7 @@
scene->forces.addTorque(id2,(geom->radius2+.5*(phys->refPD-geom->penetrationDepth))*geom->normal.cross(f));
}
TIMING_DELTAS_CHECKPOINT("rest");
+ return true;
}
=== modified file 'pkg/dem/ConcretePM.hpp'
--- pkg/dem/ConcretePM.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/ConcretePM.hpp 2014-07-18 18:18:50 +0000
@@ -255,7 +255,7 @@
class Law2_ScGeom_CpmPhys_Cpm: public LawFunctor{
public:
- void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
Real elasticEnergy();
Real yieldSigmaTMagnitude(Real sigmaN, Real omega, Real undamagedCohesion, Real tanFrictionAngle) {
=== modified file 'pkg/dem/ElasticContactLaw.cpp'
--- pkg/dem/ElasticContactLaw.cpp 2013-03-06 17:30:45 +0000
+++ pkg/dem/ElasticContactLaw.cpp 2014-07-18 18:18:50 +0000
@@ -46,7 +46,7 @@
}
CREATE_LOGGER(Law2_ScGeom_FrictPhys_CundallStrack);
-void Law2_ScGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGeom_FrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
ScGeom* geom= static_cast<ScGeom*>(ig.get());
@@ -55,8 +55,8 @@
if (neverErase) {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();}
- else scene->interactions->requestErase(contact);
- return;}
+ else return false;
+ }
Real& un=geom->penetrationDepth;
phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
@@ -98,7 +98,7 @@
}
}
-void Law2_ScGeom_ViscoFrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGeom_ViscoFrictPhys_CundallStrack::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
ScGeom* geom= static_cast<ScGeom*>(ig.get());
ViscoFrictPhys* phys = static_cast<ViscoFrictPhys*>(ip.get());
if (shearCreep) {
=== modified file 'pkg/dem/ElasticContactLaw.hpp'
--- pkg/dem/ElasticContactLaw.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/ElasticContactLaw.hpp 2014-07-18 18:18:50 +0000
@@ -16,7 +16,7 @@
class Law2_ScGeom_FrictPhys_CundallStrack: public LawFunctor{
public:
OpenMPAccumulator<Real> plasticDissipation;
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
Real elasticEnergy ();
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
@@ -38,7 +38,7 @@
class Law2_ScGeom_ViscoFrictPhys_CundallStrack: public Law2_ScGeom_FrictPhys_CundallStrack{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_ViscoFrictPhys_CundallStrack,Law2_ScGeom_FrictPhys_CundallStrack,"Law similar to :yref:`Law2_ScGeom_FrictPhys_CundallStrack` with the addition of shear creep at contacts.",
((bool,shearCreep,false,," "))
((Real,viscosity,1,," "))
=== modified file 'pkg/dem/FrictViscoPM.cpp'
--- pkg/dem/FrictViscoPM.cpp 2014-01-24 21:59:41 +0000
+++ pkg/dem/FrictViscoPM.cpp 2014-07-18 18:18:50 +0000
@@ -154,7 +154,7 @@
CREATE_LOGGER(Law2_ScGeom_FrictViscoPhys_CundallStrackVisco);
-void Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
+bool Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
{
LOG_TRACE( "Law2_ScGeom_FrictViscoPhys_CundallStrackVisco::go - create interaction physics" );
@@ -166,8 +166,7 @@
if (neverErase) {
phys->shearForce = Vector3r::Zero();
phys->normalForce = Vector3r::Zero();}
- else scene->interactions->requestErase(contact);
- return;
+ else return false;
}
Real& un=geom->penetrationDepth;
phys->normalForce = phys->kn*std::max(un,(Real) 0) * geom->normal;
@@ -224,5 +223,6 @@
scene->forces.addTorque(id1,(geom->radius1-0.5*geom->penetrationDepth)* geom->normal.cross(force));
scene->forces.addTorque(id2,(geom->radius2-0.5*geom->penetrationDepth)* geom->normal.cross(force));
}
+ return true;
}
=== modified file 'pkg/dem/FrictViscoPM.hpp'
--- pkg/dem/FrictViscoPM.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/FrictViscoPM.hpp 2014-07-18 18:18:50 +0000
@@ -95,7 +95,7 @@
class Law2_ScGeom_FrictViscoPhys_CundallStrackVisco: public LawFunctor{
public:
OpenMPAccumulator<Real> plasticDissipation;
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
Real elasticEnergy ();
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
=== modified file 'pkg/dem/HertzMindlin.cpp'
--- pkg/dem/HertzMindlin.cpp 2014-07-14 20:08:33 +0000
+++ pkg/dem/HertzMindlin.cpp 2014-07-18 18:18:50 +0000
@@ -146,20 +146,20 @@
return adhesionEnergy;
}
-void Law2_ScGeom_MindlinPhys_MindlinDeresiewitz::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGeom_MindlinPhys_MindlinDeresiewitz::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
Body::id_t id1(contact->getId1()), id2(contact->getId2());
ScGeom* geom = static_cast<ScGeom*>(ig.get());
MindlinPhys* phys=static_cast<MindlinPhys*>(ip.get());
const Real uN=geom->penetrationDepth;
if (uN<0) {
- if (neverErase) {phys->shearForce = phys->normalForce = Vector3r::Zero(); phys->kn=phys->ks=0; return;}
- else {scene->interactions->requestErase(contact); return;}
+ if (neverErase) {phys->shearForce = phys->normalForce = Vector3r::Zero(); phys->kn=phys->ks=0; return true;}
+ else {return false;}
}
// normal force
Real Fn=phys->kno*pow(uN,3/2.);
phys->normalForce=Fn*geom->normal;
// exactly zero would not work with the shear formulation, and would give zero shear force anyway
- if(Fn==0) return;
+ if(Fn==0) return true;
//phys->kn=3./2.*phys->kno*std::pow(uN,0.5); // update stiffness, not needed
// contact radius
@@ -188,14 +188,14 @@
scene->forces.addTorque(id2,(geom->radius2-.5*geom->penetrationDepth)*geom->normal.cross(f));
}
-void Law2_ScGeom_MindlinPhys_HertzWithLinearShear::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGeom_MindlinPhys_HertzWithLinearShear::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
Body::id_t id1(contact->getId1()), id2(contact->getId2());
ScGeom* geom = static_cast<ScGeom*>(ig.get());
MindlinPhys* phys=static_cast<MindlinPhys*>(ip.get());
const Real uN=geom->penetrationDepth;
if (uN<0) {
- if (neverErase) {phys->shearForce = phys->normalForce = Vector3r::Zero(); phys->kn=phys->ks=0; return;}
- else {scene->interactions->requestErase(contact); return;}
+ if (neverErase) {phys->shearForce = phys->normalForce = Vector3r::Zero(); phys->kn=phys->ks=0; return true;}
+ else return false;
}
// normal force
Real Fn=phys->kno*pow(uN,3/2.);
@@ -228,13 +228,14 @@
scene->forces.addForce(id2,-f);
scene->forces.addTorque(id1,(geom->radius1-.5*geom->penetrationDepth)*geom->normal.cross(f));
scene->forces.addTorque(id2,(geom->radius2-.5*geom->penetrationDepth)*geom->normal.cross(f));
+ return true;
}
/******************** Law2_ScGeom_MindlinPhys_Mindlin *********/
CREATE_LOGGER(Law2_ScGeom_MindlinPhys_Mindlin);
-void Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
const Real& dt = scene->dt; // get time step
Body::id_t id1 = contact->getId1(); // get id body 1
@@ -262,8 +263,8 @@
Real uN = scg->penetrationDepth; // get overlapping
if (uN<0) {
- if (neverErase) {phys->shearForce = phys->normalForce = Vector3r::Zero(); phys->kn=phys->ks=0; return;}
- else {scene->interactions->requestErase(contact); return;}
+ if (neverErase) {phys->shearForce = phys->normalForce = Vector3r::Zero(); phys->kn=phys->ks=0; return true;}
+ else return false;
}
/* Hertz-Mindlin's formulation (PFC)
Note that the normal stiffness here is a secant value (so as it is cannot be used in the GSTS)
@@ -518,7 +519,7 @@
scene->forces.addTorque(id1,-moment);
scene->forces.addTorque(id2,moment);
}
-
+ return true;
// update variables
//phys->prevNormal = scg->normal;
}
=== modified file 'pkg/dem/HertzMindlin.hpp'
--- pkg/dem/HertzMindlin.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/HertzMindlin.hpp 2014-07-18 18:18:50 +0000
@@ -93,7 +93,7 @@
class Law2_ScGeom_MindlinPhys_MindlinDeresiewitz: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
FUNCTOR2D(ScGeom,MindlinPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_MindlinPhys_MindlinDeresiewitz,LawFunctor,
"Hertz-Mindlin contact law with partial slip solution, as described in [Thornton1991]_.",
@@ -104,7 +104,7 @@
class Law2_ScGeom_MindlinPhys_HertzWithLinearShear: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
FUNCTOR2D(ScGeom,MindlinPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_MindlinPhys_HertzWithLinearShear,LawFunctor,
"Constitutive law for the Hertz formulation (using :yref:`MindlinPhys.kno`) and linear beahvior in shear (using :yref:`MindlinPhys.kso` for stiffness and :yref:`FrictPhys.tangensOfFrictionAngle`). \n\n.. note:: No viscosity or damping. If you need those, look at :yref:`Law2_ScGeom_MindlinPhys_Mindlin`, which also includes non-linear Mindlin shear.",
@@ -119,7 +119,7 @@
class Law2_ScGeom_MindlinPhys_Mindlin: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
Real normElastEnergy();
Real adhesionEnergy();
=== modified file 'pkg/dem/InelastCohFrictPM.cpp'
--- pkg/dem/InelastCohFrictPM.cpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/InelastCohFrictPM.cpp 2014-07-18 18:18:50 +0000
@@ -91,7 +91,7 @@
}
-void Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
+bool Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact)
{
//FIXME : non cohesive contact are not implemented, it would be useful to use setCohesionNow, setCohesionOnNewContacts etc ...
const int &id1 = contact->getId1();
@@ -112,8 +112,7 @@
if(-un>phys->maxExten || phys->isBroken){//plastic failure.
phys->isBroken=1;
phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
- scene->interactions->requestErase(contact);
- return;
+ return false;
}
Fn=phys->knT*un; //elasticity
if(-Fn>phys->maxElT || phys->onPlastT){ //so we are on plastic deformation.
@@ -144,9 +143,9 @@
phys->isBroken=1;
phys->normalForce=phys->shearForce=phys->moment_twist=phys->moment_bending=Vector3r(0,0,0);
if(geom->penetrationDepth<=0){ //do not erase the contact while penetrationDepth<0 because it would be recreated at next timestep.
- scene->interactions->requestErase(contact);
+ return false;
}
- return;
+ return true;
}
Fn=phys->knC*un;
if(Fn>phys->maxElC || phys->onPlastC){
@@ -268,4 +267,5 @@
applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position + (scene->isPeriodic ? scene->cell->intrShiftPos(contact->cellDist): Vector3r::Zero()));
scene->forces.addTorque(id1,-phys->moment_bending-phys->moment_twist);
scene->forces.addTorque(id2,phys->moment_bending+phys->moment_twist);
+ return true;
}
=== modified file 'pkg/dem/InelastCohFrictPM.hpp'
--- pkg/dem/InelastCohFrictPM.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/InelastCohFrictPM.hpp 2014-07-18 18:18:50 +0000
@@ -137,7 +137,7 @@
public:
Real normElastEnergy();
Real shearElastEnergy();
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment,LawFunctor,"This law is currently under developpement. Final version and documentation will come before the end of 2014.",
,,
.def("normElastEnergy",&Law2_ScGeom6D_InelastCohFrictPhys_CohesionMoment::normElastEnergy,"Compute normal elastic energy.")
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.cpp'
--- pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.cpp 2014-07-18 18:18:50 +0000
@@ -11,7 +11,7 @@
/********************** Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM ****************************/
CREATE_LOGGER(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM);
-void Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
const int &id1 = contact->getId1();
const int &id2 = contact->getId2();
@@ -42,7 +42,7 @@
if ( smoothJoint && phys->isOnJoint ) {
if ( phys->more || ( phys->jointCumulativeSliding > (2*min(geom->radius1,geom->radius2)) ) ) {
- scene->interactions->requestErase(contact); return;
+ return false;
} else {
D = phys->initD - std::abs((b1->state->pos - b2->state->pos).dot(phys->jointNormal));
}
@@ -54,7 +54,7 @@
if (D < 0) { //spheres do not "touch" !
if ( !phys->isCohesive )
{
- scene->interactions->requestErase(contact); return;
+ return false;
}
if ( phys->isCohesive && (phys->FnMax>0) && (std::abs(D)>Dtensile) ) {
@@ -78,7 +78,7 @@
cracksFileExist=true;
// delete contact
- scene->interactions->requestErase(contact); return;
+ return false;
}
}
@@ -143,7 +143,7 @@
// delete contact if in tension, set the contact properties to friction if in compression
if ( D < 0 )
{
- scene->interactions->requestErase(contact); return;
+ return false;
}
else
{
@@ -165,11 +165,12 @@
scene->forces.addForce (id2, f);
// simple solution to avoid torque computation for particles interacting on a smooth joint
- if ( (phys->isOnJoint)&&(smoothJoint) ) return;
+ if ( (phys->isOnJoint)&&(smoothJoint) ) return true;
/// those lines are needed if rootBody->forces.addForce and rootBody->forces.addMoment are used instead of applyForceAtContactPoint -> NOTE need to check for accuracy!!!
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));
+ return true;
}
=== modified file 'pkg/dem/JointedCohesiveFrictionalPM.hpp'
--- pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-06-10 14:27:43 +0000
+++ pkg/dem/JointedCohesiveFrictionalPM.hpp 2014-07-18 18:18:50 +0000
@@ -95,7 +95,7 @@
/** 2d functor creating the interaction law (Law2) based on SphereContactGeometry (ScGeom) and JCFpmPhys of 2 bodies, returning type JointedCohesiveFrictionalPM */
class Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
FUNCTOR2D(ScGeom,JCFpmPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM,LawFunctor,"Interaction law for cohesive frictional material, e.g. rock, possibly presenting joint surfaces, that can be mechanically described with a smooth contact logic [Ivars2011]_ (implemented in Yade in [Scholtes2012]_). See examples/jointedCohesiveFrictionalPM for script examples. Joint surface definitions (through stl meshes or direct definition with gts module) are illustrated there.",
=== modified file 'pkg/dem/L3Geom.cpp'
--- pkg/dem/L3Geom.cpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/L3Geom.cpp 2014-07-18 18:18:50 +0000
@@ -284,14 +284,14 @@
return true;
}
-void Law2_L3Geom_FrictPhys_ElPerfPl::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
+bool Law2_L3Geom_FrictPhys_ElPerfPl::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
L3Geom* geom=static_cast<L3Geom*>(ig.get()); FrictPhys* phys=static_cast<FrictPhys*>(ip.get());
// compute force
Vector3r& localF(geom->F);
localF=geom->relU().cwiseProduct(Vector3r(phys->kn,phys->ks,phys->ks));
// break if necessary
- if(localF[0]>0 && !noBreak){ scene->interactions->requestErase(I); return; }
+ if(localF[0]>0 && !noBreak) return false;
if(!noSlip){
// plastic slip, if necessary; non-zero elastic limit only for compression
@@ -308,10 +308,11 @@
if(scene->trackEnergy) { scene->energy->add(0.5*(pow(geom->relU()[0],2)*phys->kn+(pow(geom->relU()[1],2)+pow(geom->relU()[2],2))*phys->ks),"elastPotential",elastPotentialIx,/*reset at every timestep*/true); }
// apply force: this converts the force to global space, updates NormShearPhys::{normal,shear}Force, applies to particles
geom->applyLocalForce(localF,I,scene,phys);
+ return true;
}
-void Law2_L6Geom_FrictPhys_Linear::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
+bool Law2_L6Geom_FrictPhys_Linear::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
L6Geom& geom=ig->cast<L6Geom>(); FrictPhys& phys=ip->cast<FrictPhys>();
// simple linear relationships
=== modified file 'pkg/dem/L3Geom.hpp'
--- pkg/dem/L3Geom.hpp 2014-01-20 16:14:58 +0000
+++ pkg/dem/L3Geom.hpp 2014-07-18 18:18:50 +0000
@@ -174,7 +174,7 @@
struct Law2_L3Geom_FrictPhys_ElPerfPl: public LawFunctor{
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
FUNCTOR2D(L3Geom,FrictPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_L3Geom_FrictPhys_ElPerfPl,LawFunctor,"Basic law for testing :yref:`L3Geom`; it bears no cohesion (unless *noBreak* is ``True``), and plastic slip obeys the Mohr-Coulomb criterion (unless *noSlip* is ``True``).",
((bool,noBreak,false,,"Do not break contacts when particles separate."))
@@ -186,7 +186,7 @@
REGISTER_SERIALIZABLE(Law2_L3Geom_FrictPhys_ElPerfPl);
struct Law2_L6Geom_FrictPhys_Linear: public Law2_L3Geom_FrictPhys_ElPerfPl{
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
FUNCTOR2D(L6Geom,FrictPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_L6Geom_FrictPhys_Linear,Law2_L3Geom_FrictPhys_ElPerfPl,"Basic law for testing :yref:`L6Geom` -- linear in both normal and shear sense, without slip or breakage.",
((Real,charLen,1,,"Characteristic length with the meaning of the stiffness ratios bending/shear and torsion/normal."))
=== modified file 'pkg/dem/LudingPM.cpp'
--- pkg/dem/LudingPM.cpp 2013-10-11 19:22:48 +0000
+++ pkg/dem/LudingPM.cpp 2014-07-18 18:18:50 +0000
@@ -71,7 +71,7 @@
return 2.0*a;
}
-void Law2_ScGeom_LudingPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
+bool Law2_ScGeom_LudingPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I){
const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
LudingPhys& phys=*static_cast<LudingPhys*>(_phys.get());
@@ -82,8 +82,7 @@
const Real Delt = geom.penetrationDepth;
if (Delt<0) {
- scene->interactions->requestErase(I);
- return;
+ return false;
};
const BodyContainer& bodies = *scene->bodies;
@@ -195,6 +194,6 @@
addTorque(id1,-c1x.cross(f),scene);
addTorque(id2, c2x.cross(f),scene);
}
-
+ return true;
}
=== modified file 'pkg/dem/LudingPM.hpp'
--- pkg/dem/LudingPM.hpp 2013-10-11 14:42:51 +0000
+++ pkg/dem/LudingPM.hpp 2014-07-18 18:18:50 +0000
@@ -59,7 +59,7 @@
class Law2_ScGeom_LudingPhys_Basic: public LawFunctor {
public :
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
private:
Real calculateCapillarForce(const ScGeom& geom, LudingPhys& phys);
FUNCTOR2D(ScGeom,LudingPhys);
=== modified file 'pkg/dem/NormalInelasticPM.cpp'
--- pkg/dem/NormalInelasticPM.cpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/NormalInelasticPM.cpp 2014-07-18 18:18:50 +0000
@@ -14,7 +14,7 @@
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)
+bool Law2_ScGeom6D_NormalInelasticityPhys_NormalInelasticity::go(shared_ptr<IGeom>& iG, shared_ptr<IPhys>& iP, Interaction* contact)
{
int id1 = contact->getId1();
int id2 = contact->getId2();
@@ -39,8 +39,7 @@
// 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;
+ return false;// this, among other things, resets the interaction : geometry and physics variables (as forces, ...) are reset to defaut values
}
@@ -87,10 +86,7 @@
// ******** 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();
+ return false;
}
else
{
@@ -145,6 +141,7 @@
}
// ******** Moment law END ******* //
}
+ return true;
}
=== modified file 'pkg/dem/NormalInelasticPM.hpp'
--- pkg/dem/NormalInelasticPM.hpp 2014-07-14 20:08:34 +0000
+++ pkg/dem/NormalInelasticPM.hpp 2014-07-18 18:18:50 +0000
@@ -66,7 +66,7 @@
,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*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
FUNCTOR2D(ScGeom,NormalInelasticityPhys);
=== modified file 'pkg/dem/Polyhedra.cpp'
--- pkg/dem/Polyhedra.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/Polyhedra.cpp 2014-07-18 18:18:50 +0000
@@ -413,11 +413,11 @@
//**************************************************************************************
// Apply forces on polyhedrons in collision based on geometric configuration
-void PolyhedraVolumetricLaw::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
+bool PolyhedraVolumetricLaw::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* I){
- if (!I->geom) {return;}
+ if (!I->geom) {return true;}
const shared_ptr<PolyhedraGeom>& contactGeom(YADE_PTR_DYN_CAST<PolyhedraGeom>(I->geom));
- if(!contactGeom) {return;}
+ if(!contactGeom) {return true;}
const Body::id_t idA=I->getId1(), idB=I->getId2();
const shared_ptr<Body>& A=Body::byId(idA), B=Body::byId(idB);
@@ -425,13 +425,11 @@
//erase the interaction when aAbB shows separation, otherwise keep it to be able to store previous separating plane for fast detection of separation
if (A->bound->min[0] >= B->bound->max[0] || B->bound->min[0] >= A->bound->max[0] || A->bound->min[1] >= B->bound->max[1] || B->bound->min[1] >= A->bound->max[1] || A->bound->min[2] >= B->bound->max[2] || B->bound->min[2] >= A->bound->max[2]) {
- scene->interactions->requestErase(I);
- phys->normalForce = Vector3r(0.,0.,0.); phys->shearForce = Vector3r(0.,0.,0.);
- return;
+ return false;
}
//zero penetration depth means no interaction force
- if(!(contactGeom->equivalentPenetrationDepth > 1E-18) || !(contactGeom->penetrationVolume > 0)) {phys->normalForce = Vector3r(0.,0.,0.); phys->shearForce = Vector3r(0.,0.,0.); return;}
+ if(!(contactGeom->equivalentPenetrationDepth > 1E-18) || !(contactGeom->penetrationVolume > 0)) {phys->normalForce = Vector3r(0.,0.,0.); phys->shearForce = Vector3r(0.,0.,0.); return true;}
Vector3r normalForce=contactGeom->normal*contactGeom->penetrationVolume*phys->kn;
//shear force: in case the polyhdras are separated and come to contact again, one should not use the previous shear force
@@ -491,6 +489,7 @@
//needed to be able to acces interaction forces in other parts of yade
phys->normalForce = normalForce;
phys->shearForce = shearForce;
+ return true;
}
#endif // YADE_CGAL
=== modified file 'pkg/dem/Polyhedra.hpp'
--- pkg/dem/Polyhedra.hpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/Polyhedra.hpp 2014-07-18 18:18:50 +0000
@@ -242,7 +242,7 @@
class PolyhedraVolumetricLaw: public LawFunctor{
OpenMPAccumulator<Real> plasticDissipation;
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
Real elasticEnergy ();
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
=== modified file 'pkg/dem/Tetra.cpp'
--- pkg/dem/Tetra.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/Tetra.cpp 2014-07-18 18:18:50 +0000
@@ -1208,13 +1208,12 @@
#ifdef YADE_CGAL
-void Law2_TTetraSimpleGeom_NormPhys_Simple::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_TTetraSimpleGeom_NormPhys_Simple::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
int id1 = contact->getId1(), id2 = contact->getId2();
TTetraSimpleGeom* geom= static_cast<TTetraSimpleGeom*>(ig.get());
NormPhys* phys = static_cast<NormPhys*>(ip.get());
if ( geom->flag == 0 || geom->penetrationVolume <= 0. ) {
- scene->interactions->requestErase(contact);
- return;
+ return false;
}
Real& un=geom->penetrationVolume;
phys->normalForce=phys->kn*std::max(un,(Real) 0)*geom->normal;
@@ -1223,6 +1222,7 @@
State* de2 = Body::byId(id2,scene)->state.get();
applyForceAtContactPoint(-phys->normalForce, geom->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
// TODO periodic
+ return true;
}
#endif
=== modified file 'pkg/dem/Tetra.hpp'
--- pkg/dem/Tetra.hpp 2014-07-03 17:20:40 +0000
+++ pkg/dem/Tetra.hpp 2014-07-18 18:18:50 +0000
@@ -169,7 +169,7 @@
class Law2_TTetraSimpleGeom_NormPhys_Simple: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC(Law2_TTetraSimpleGeom_NormPhys_Simple,LawFunctor,"EXPERIMENTAL. TODO");
FUNCTOR2D(TTetraSimpleGeom,NormPhys);
DECLARE_LOGGER;
=== modified file 'pkg/dem/ViscoelasticCapillarPM.cpp'
--- pkg/dem/ViscoelasticCapillarPM.cpp 2014-06-18 07:35:57 +0000
+++ pkg/dem/ViscoelasticCapillarPM.cpp 2014-07-18 18:18:50 +0000
@@ -62,7 +62,7 @@
}
/* Law2_ScGeom_ViscElCapPhys_Basic */
-void Law2_ScGeom_ViscElCapPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I) {
+bool Law2_ScGeom_ViscElCapPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I) {
Vector3r force = Vector3r::Zero();
const int id1 = I->getId1();
@@ -111,7 +111,7 @@
addForce (id1,-phys.normalForce,scene);
addForce (id2, phys.normalForce,scene);
};
- return;
+ return true;
} else {
if (phys.liqBridgeActive) {
VLiqBridg -= phys.Vb;
@@ -127,8 +127,7 @@
scene->delIntrs.push_back(B1);
scene->delIntrs.push_back(B2);
#endif
- scene->interactions->requestErase(I);
- return;
+ return false;
};
};
@@ -150,6 +149,7 @@
addTorque(id1, torque1,scene);
addTorque(id2, torque2,scene);
}
+ return true;
}
Real Law2_ScGeom_ViscElCapPhys_Basic::critDist(const Real& Vb, const Real& R, const Real& Theta) {
=== modified file 'pkg/dem/ViscoelasticCapillarPM.hpp'
--- pkg/dem/ViscoelasticCapillarPM.hpp 2014-06-25 14:46:14 +0000
+++ pkg/dem/ViscoelasticCapillarPM.hpp 2014-07-18 18:18:50 +0000
@@ -61,7 +61,7 @@
/// Constitutive law
class Law2_ScGeom_ViscElCapPhys_Basic: public LawFunctor {
public :
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
static Real Willett_numeric_f (const ScGeom& geom, ViscElCapPhys& phys);
static Real Willett_analytic_f (const ScGeom& geom, ViscElCapPhys& phys);
static Real Weigert_f (const ScGeom& geom, ViscElCapPhys& phys);
=== modified file 'pkg/dem/ViscoelasticPM.cpp'
--- pkg/dem/ViscoelasticPM.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/ViscoelasticPM.cpp 2014-07-18 18:18:50 +0000
@@ -29,7 +29,7 @@
}
/* Law2_ScGeom_ViscElPhys_Basic */
-void Law2_ScGeom_ViscElPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I) {
+bool Law2_ScGeom_ViscElPhys_Basic::go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I) {
Vector3r force = Vector3r::Zero();
Vector3r torque1 = Vector3r::Zero();
Vector3r torque2 = Vector3r::Zero();
@@ -41,11 +41,8 @@
addForce (id2, force,scene);
addTorque(id1, torque1,scene);
addTorque(id2, torque2,scene);
- return;
- } else {
- scene->interactions->requestErase(I);
- return;
- }
+ return true;
+ } else return false;
}
bool computeForceTorqueViscEl(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I, Vector3r & force, Vector3r & torque1, Vector3r & torque2) {
=== modified file 'pkg/dem/ViscoelasticPM.hpp'
--- pkg/dem/ViscoelasticPM.hpp 2014-06-03 08:53:31 +0000
+++ pkg/dem/ViscoelasticPM.hpp 2014-07-18 18:18:50 +0000
@@ -90,7 +90,7 @@
/// This class provides linear viscoelastic contact model
class Law2_ScGeom_ViscElPhys_Basic: public LawFunctor {
public :
- virtual void go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
+ virtual bool go(shared_ptr<IGeom>&, shared_ptr<IPhys>&, Interaction*);
public :
FUNCTOR2D(ScGeom,ViscElPhys);
YADE_CLASS_BASE_DOC(Law2_ScGeom_ViscElPhys_Basic,LawFunctor,"Linear viscoelastic model operating on :yref:`ScGeom` and :yref:`ViscElPhys`. The model is mostly based on the paper for For details see Pournin [Pournin2001]_ .");
=== modified file 'pkg/dem/WirePM.cpp'
--- pkg/dem/WirePM.cpp 2014-07-03 07:16:58 +0000
+++ pkg/dem/WirePM.cpp 2014-07-18 18:18:50 +0000
@@ -80,7 +80,7 @@
/********************** Law2_ScGeom_WirePhys_WirePM ****************************/
CREATE_LOGGER(Law2_ScGeom_WirePhys_WirePM);
-void Law2_ScGeom_WirePhys_WirePM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
+bool Law2_ScGeom_WirePhys_WirePM::go(shared_ptr<IGeom>& ig, shared_ptr<IPhys>& ip, Interaction* contact){
LOG_TRACE( "Law2_ScGeom_WirePhys_WirePM::go - contact law" );
@@ -102,8 +102,7 @@
/* check whether the particles are linked or not */
if ( !phys->isLinked ) { // destroy the interaction before calculation
- scene->interactions->requestErase(contact);
- return;
+ return false;
}
if ( (phys->isLinked) && (D < DFValues.back()(0)) ) { // spheres are linked but failure because of reaching maximal admissible displacement
phys->isLinked=false;
@@ -112,8 +111,7 @@
WireState* st2=dynamic_cast<WireState*>(b2->state.get());
st1->numBrokenLinks+=1;
st2->numBrokenLinks+=1;
- scene->interactions->requestErase(contact);
- return;
+ return false;
}
/* compute normal force Fn */
@@ -163,6 +161,7 @@
/* set shear force to zero */
phys->shearForce = Vector3r::Zero();
+ return true;
}
=== modified file 'pkg/dem/WirePM.hpp'
--- pkg/dem/WirePM.hpp 2013-08-14 08:01:55 +0000
+++ pkg/dem/WirePM.hpp 2014-07-18 18:18:50 +0000
@@ -116,7 +116,7 @@
/** 2d functor creating the interaction law (Law2) based on SphereContactGeometry (ScGeom) and WirePhys of 2 bodies, returning type WirePM */
class Law2_ScGeom_WirePhys_WirePM: public LawFunctor{
public:
- virtual void go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
+ virtual bool go(shared_ptr<IGeom>& _geom, shared_ptr<IPhys>& _phys, Interaction* I);
FUNCTOR2D(ScGeom,WirePhys);
