yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05430
[Branch ~yade-dev/yade/trunk] Rev 2374: 1. Do not uselessly pass Scene* pointer to Law2 functors (already in LawFunctor)
------------------------------------------------------------
revno: 2374
committer: Václav Šmilauer <eudoxos@xxxxxxxx>
branch nick: trunk
timestamp: Sat 2010-07-17 19:37:54 +0200
message:
1. Do not uselessly pass Scene* pointer to Law2 functors (already in LawFunctor)
2. Attempt at Law2_ScGeom_CpmPhys_Cpm (not tested)
modified:
pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp
pkg/common/Engine/Dispatcher/LawDispatcher.cpp
pkg/common/Engine/Dispatcher/LawDispatcher.hpp
pkg/common/Engine/Dispatcher/LawFunctor.hpp
pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp
pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp
pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp
pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp
pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp
pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.hpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp
pkg/dem/Engine/GlobalEngine/CundallStrack.cpp
pkg/dem/Engine/GlobalEngine/CundallStrack.hpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp
pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp
pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp
pkg/dem/meta/ConcretePM.cpp
pkg/dem/meta/ConcretePM.hpp
pkg/dem/meta/RockPM.cpp
pkg/dem/meta/RockPM.hpp
pkg/dem/meta/Shop.cpp
pkg/dem/meta/Shop.hpp
pkg/dem/meta/ViscoelasticPM.cpp
pkg/dem/meta/ViscoelasticPM.hpp
py/plot.py
--
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/common/Engine/Dispatcher/InteractionDispatchers.cpp'
--- pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2010-07-06 11:47:43 +0000
+++ pkg/common/Engine/Dispatcher/InteractionDispatchers.cpp 2010-07-17 17:37:54 +0000
@@ -150,7 +150,7 @@
assert(!swap); // reverse call would make no sense, as the arguments are of different types
}
assert(I->functorCache.constLaw);
- I->functorCache.constLaw->go(I->interactionGeometry,I->interactionPhysics,I.get(),scene);
+ I->functorCache.constLaw->go(I->interactionGeometry,I->interactionPhysics,I.get());
// process callbacks for this interaction
if(!I->isReal()) continue; // it is possible that Law2_ functor called requestErase, hence this check
=== modified file 'pkg/common/Engine/Dispatcher/LawDispatcher.cpp'
--- pkg/common/Engine/Dispatcher/LawDispatcher.cpp 2010-03-20 12:40:44 +0000
+++ pkg/common/Engine/Dispatcher/LawDispatcher.cpp 2010-07-17 17:37:54 +0000
@@ -14,7 +14,7 @@
#endif
if(I->isReal()){
assert(I->interactionGeometry); assert(I->interactionPhysics);
- operator()(I->interactionGeometry,I->interactionPhysics,I.get(),scene);
+ operator()(I->interactionGeometry,I->interactionPhysics,I.get());
if(!I->isReal() && I->isFresh(scene)) LOG_ERROR("Law functor deleted interaction that was just created. Please report bug: either this message is spurious, or the functor (or something else) is buggy.");
}
}
=== modified file 'pkg/common/Engine/Dispatcher/LawDispatcher.hpp'
--- pkg/common/Engine/Dispatcher/LawDispatcher.hpp 2010-03-20 12:40:44 +0000
+++ pkg/common/Engine/Dispatcher/LawDispatcher.hpp 2010-07-17 17:37:54 +0000
@@ -11,7 +11,7 @@
InteractionPhysics, // 2nd base classe for dispatch
LawFunctor, // functor base class
void, // return type
- TYPELIST_4(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, Scene*),
+ TYPELIST_3(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*),
false // autosymmetry
>{
public:
=== modified file 'pkg/common/Engine/Dispatcher/LawFunctor.hpp'
--- pkg/common/Engine/Dispatcher/LawFunctor.hpp 2010-04-25 15:46:26 +0000
+++ pkg/common/Engine/Dispatcher/LawFunctor.hpp 2010-07-17 17:37:54 +0000
@@ -5,7 +5,7 @@
#include<yade/core/Scene.hpp>
class LawFunctor: public Functor2D <
- void, TYPELIST_4(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, Scene*)
+ void, TYPELIST_3(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*)
>{
public:
virtual ~LawFunctor();
@@ -13,11 +13,11 @@
void addForce (const body_id_t id, const Vector3r& f,Scene* rb){rb->forces.addForce (id,f);}
void addTorque(const body_id_t id, const Vector3r& t,Scene* rb){rb->forces.addTorque(id,t);}
/*! Convenience function to apply force and torque from one force at contact point. Not sure if this is the right place for it. */
- void applyForceAtContactPoint(const Vector3r& force, const Vector3r& contactPoint, const body_id_t id1, const Vector3r& pos1, const body_id_t id2, const Vector3r& pos2, Scene* rb){
- addForce(id1,force,rb);
- addForce(id2,-force,rb);
- addTorque(id1,(contactPoint-pos1).cross(force),rb);
- addTorque(id2,-(contactPoint-pos2).cross(force),rb);
+ void applyForceAtContactPoint(const Vector3r& force, const Vector3r& contactPoint, const body_id_t id1, const Vector3r& pos1, const body_id_t id2, const Vector3r& pos2){
+ addForce(id1,force,scene);
+ addForce(id2,-force,scene);
+ addTorque(id1,(contactPoint-pos1).cross(force),scene);
+ addTorque(id2,-(contactPoint-pos2).cross(force),scene);
}
YADE_CLASS_BASE_DOC(LawFunctor,Functor,"Functor for applying constitutive laws on :yref:`interactions<Interaction>`.");
};
=== modified file 'pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp'
--- pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp 2010-01-10 09:09:32 +0000
+++ pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.cpp 2010-07-17 17:37:54 +0000
@@ -3,7 +3,7 @@
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/pkg-common/NormShearPhys.hpp>
YADE_PLUGIN((ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw));
-void ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody){
+void ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I){
Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(_geom.get());
NormShearPhys* phys=static_cast<NormShearPhys*>(_phys.get());
@@ -11,7 +11,7 @@
phys->normalForce=(geom->displacementN()*phys->kn)*geom->normal;
phys->shearForce=geom->displacementT()*phys->ks;
- applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, I->getId1(), geom->se31.position, I->getId2(), geom->se32.position, rootBody);
+ applyForceAtContactPoint(phys->normalForce+phys->shearForce, geom->contactPoint, I->getId1(), geom->se31.position, I->getId2(), geom->se32.position);
}
=== modified file 'pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp'
--- pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp 2010-01-10 09:09:32 +0000
+++ pkg/dem/Engine/Functor/ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw.hpp 2010-07-17 17:37:54 +0000
@@ -3,7 +3,7 @@
/* Experimental constitutive law using the LawDispatcher.
* Has only purely elastic normal and shear components. */
class ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw: public LawFunctor {
- virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, Scene*);
+ virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*);
FUNCTOR2D(Dem3DofGeom,NormShearPhys);
REGISTER_CLASS_AND_BASE(ef2_Dem3Dof_NormalShear_ElasticFrictionalLaw,LawFunctor);
};
=== modified file 'pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp'
--- pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2010-05-13 20:19:38 +0000
+++ pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.cpp 2010-07-17 17:37:54 +0000
@@ -11,7 +11,7 @@
#include<yade/pkg-dem/ViscoelasticInteraction.hpp>
#include<yade/pkg-common/RigidBodyParameters.hpp>
YADE_PLUGIN((ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw));
-void ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody){
+void ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I){
ScGeom* geom=static_cast<ScGeom*>(_geom.get());
ViscoelasticInteraction* phys=static_cast<ViscoelasticInteraction*>(_phys.get());
@@ -20,18 +20,18 @@
int id2 = I->getId2();
if (geom->penetrationDepth<0) {
- rootBody->interactions->requestErase(id1,id2);
+ scene->interactions->requestErase(id1,id2);
return;
}
- shared_ptr<BodyContainer>& bodies = rootBody->bodies;
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
RigidBodyParameters* de1 = YADE_CAST<RigidBodyParameters*>((*bodies)[id1]->physicalParameters.get());
RigidBodyParameters* de2 = YADE_CAST<RigidBodyParameters*>((*bodies)[id2]->physicalParameters.get());
Vector3r& shearForce = phys->shearForce;
- if (I->isFresh(rootBody)) shearForce=Vector3r(0,0,0);
+ if (I->isFresh(scene)) shearForce=Vector3r(0,0,0);
const Real& dt = scene->dt;
@@ -67,10 +67,10 @@
}
Vector3r f = phys->normalForce + shearForce;
- addForce (id1,-f,rootBody);
- addForce (id2, f,rootBody);
- addTorque(id1,-c1x.Cross(f),rootBody);
- addTorque(id2, c2x.Cross(f),rootBody);
+ addForce (id1,-f,scene);
+ addForce (id2, f,scene);
+ addTorque(id1,-c1x.Cross(f),scene);
+ addTorque(id2, c2x.Cross(f),scene);
}
YADE_REQUIRE_FEATURE(PHYSPAR);
=== modified file 'pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp'
--- pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp 2009-12-13 20:11:31 +0000
+++ pkg/dem/Engine/Functor/ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw.hpp 2010-07-17 17:37:54 +0000
@@ -14,7 +14,7 @@
class ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw : public LawFunctor
{
public :
- virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, Scene*);
+ virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*);
FUNCTOR2D(ScGeom,ViscoelasticInteraction);
REGISTER_CLASS_AND_BASE(ef2_Spheres_Viscoelastic_SimpleViscoelasticContactLaw,LawFunctor);
};
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 2010-07-12 09:11:16 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.cpp 2010-07-17 17:37:54 +0000
@@ -16,13 +16,13 @@
/********************** Law2_SCG_MomentPhys_CohesionlessMomentRotation ****************************/
CREATE_LOGGER(Law2_SCG_MomentPhys_CohesionlessMomentRotation);
-void Law2_SCG_MomentPhys_CohesionlessMomentRotation::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* rootBody){
+void Law2_SCG_MomentPhys_CohesionlessMomentRotation::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
ScGeom* geom = static_cast<ScGeom*>(ig.get()); //InteractionGeometry
MomentPhys* phys = static_cast<MomentPhys*>(ip.get()); //InteractionPhysics
int id1 = contact->getId1(); //Id of body1
int id2 = contact->getId2(); //Id of body2
- shared_ptr<BodyContainer>& bodies = rootBody->bodies;
+ shared_ptr<BodyContainer>& bodies = scene->bodies;
Body* b1 = ( *bodies ) [id1].get();
Body* b2 = ( *bodies ) [id2].get();
const Real& dt = scene->dt; //Get TimeStep
@@ -30,7 +30,7 @@
/*NormalForce */
Real displN=geom->penetrationDepth; //get penetrationDepth from SCG
- if (displN<0){rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;} //! NOTE: the sign for penetrationdepth is different from SCG and Dem3DofGeom
+ if (displN<0){scene->interactions->requestErase(contact->getId1(),contact->getId2()); return;} //! NOTE: the sign for penetrationdepth is different from SCG and Dem3DofGeom
Real Fn =phys->kn*displN; //scalar force
phys->normalForce= Fn*geom->normal; //vector force NOTE: normal is position2-position1. Therefore it is directed from particle1 to particle2
@@ -107,10 +107,10 @@
/* Apply forces */
Vector3r f = phys->normalForce + shearForce;
- rootBody->forces.addForce (id1,-f);
- rootBody->forces.addForce (id2, f);
- rootBody->forces.addTorque(id1,-(c1x).cross(f)); //about axis perpendicular to normal and force
- rootBody->forces.addTorque(id2, c2x.cross(f));
+ scene->forces.addForce (id1,-f);
+ scene->forces.addForce (id2, f);
+ scene->forces.addTorque(id1,-(c1x).cross(f)); //about axis perpendicular to normal and force
+ scene->forces.addTorque(id2, c2x.cross(f));
/* Moment Rotation Law */
Quaternionr delta( b1->state->ori * phys->initialOrientation1.conjugate() *phys->initialOrientation2 * b2->state->ori.conjugate()); //relative orientation
@@ -148,8 +148,8 @@
phys->moment_twist = moment_twist;
phys->moment_bending = moment_bending;
- rootBody->forces.addTorque(id1,-moment);
- rootBody->forces.addTorque(id2, moment);
+ scene->forces.addTorque(id1,-moment);
+ scene->forces.addTorque(id2, moment);
/// Moment law END ///
phys->prevNormal = geom->normal;
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.hpp'
--- pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.hpp 2010-04-25 13:18:11 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesionlessMomentRotation.hpp 2010-07-17 17:37:54 +0000
@@ -12,7 +12,7 @@
class Law2_SCG_MomentPhys_CohesionlessMomentRotation: public LawFunctor{
public:
- virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
FUNCTOR2D(ScGeom,MomentPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_SCG_MomentPhys_CohesionlessMomentRotation,LawFunctor,"Contact law based on Plassiard et al. (2009) : A spherical discrete element model: calibration procedure and incremental response. The functionality has been verified with results in the paper.\n\nThe contribution of stiffnesses are scaled according to the radius of the particle, as implemented in that paper.\n\nSee also associated classes :yref:`MomentMat`, :yref:`Ip2_MomentMat_MomentMat_MomentPhys`, :yref:`MomentPhys`.\n\n.. note::\n\tThis constitutive law can be used with triaxial test, but the following significant changes in code have to be made: :yref:`Ip2_MomentMat_MomentMat_MomentPhys` and :yref:`Law2_SCG_MomentPhys_CohesionlessMomentRotation` have to be added. Since it uses :yref:`ScGeom`, it uses :yref:`boxes<Box>` rather than :yref:`facets<Facet>`. :yref:`Spheres<Sphere>` and :yref:`boxes<Box>` have to be changed to :yref:`MomentMat` rather than :yref:`FrictMat`.",
((bool,preventGranularRatcheting,false,"??"))
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-17 17:37:54 +0000
@@ -34,12 +34,12 @@
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
if(!I->isReal()) continue;
- functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), scene);
+ functor->go(I->interactionGeometry, I->interactionPhysics, I.get());
}
}
-void Law2_ScGeom_CohFrictPhys_ElasticPlastic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb)
+void Law2_ScGeom_CohFrictPhys_ElasticPlastic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact)
{
const Real& dt = scene->dt;
// if (detectBrokenBodies //Experimental, has no effect
@@ -49,28 +49,28 @@
// YADE_CAST<CohFrictMat*> ((*bodies)[contact->getId2()]->material.get())->isBroken = false;}
const int &id1 = contact->getId1();
const int &id2 = contact->getId2();
- Body* b1 = Body::byId(id1,ncb).get();
- Body* b2 = Body::byId(id2,ncb).get();
+ Body* b1 = Body::byId(id1,scene).get();
+ Body* b2 = Body::byId(id2,scene).get();
ScGeom* currentContactGeometry = YADE_CAST<ScGeom*> (ig.get());
CohFrictPhys* currentContactPhysics = YADE_CAST<CohFrictPhys*> (ip.get());
Vector3r& shearForce = currentContactPhysics->shearForce;
- if (contact->isFresh(ncb)) shearForce = Vector3r::Zero();
+ if (contact->isFresh(scene)) shearForce = Vector3r::Zero();
Real un = currentContactGeometry->penetrationDepth;
Real Fn = currentContactPhysics->kn*un;
currentContactPhysics->normalForce = Fn*currentContactGeometry->normal;
if (un < 0 && (currentContactPhysics->normalForce.squaredNorm() > pow(currentContactPhysics->normalAdhesion,2)
|| currentContactPhysics->normalAdhesion==0)) {
// BREAK due to tension
- ncb->interactions->requestErase(contact->getId1(),contact->getId2());
+ scene->interactions->requestErase(contact->getId1(),contact->getId2());
// contact->interactionPhysics was reset now; currentContactPhysics still hold the object, but is not associated with the interaction anymore
// currentContactPhysics->cohesionBroken = true;
// currentContactPhysics->normalForce = Vector3r::ZERO;
// currentContactPhysics->shearForce = Vector3r::ZERO;
} else {
- State* de1 = Body::byId(id1,ncb)->state.get();
- State* de2 = Body::byId(id2,ncb)->state.get();
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
///////////////////////// CREEP START ///////////
if (shear_creep) shearForce -= currentContactPhysics->ks*(shearForce*dt/creep_viscosity);
///////////////////////// CREEP END ////////////
@@ -99,7 +99,7 @@
shearForce *= maxFs;
if (Fn<0) currentContactPhysics->normalForce = Vector3r::Zero();//Vector3r::Zero()
}
- applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
+ applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
/// Moment law ///
if (momentRotationLaw && (!currentContactPhysics->cohesionBroken || always_use_moment_law)) {
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp 2010-07-12 10:37:20 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.hpp 2010-07-17 17:37:54 +0000
@@ -15,7 +15,7 @@
class Law2_ScGeom_CohFrictPhys_ElasticPlastic: public LawFunctor{
public:
- virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CohFrictPhys_ElasticPlastic,LawFunctor,"Law for linear traction-compression-bending-twisting, with cohesion+friction and Mohr-Coulomb plasticity surface. Can be elastic-fragile or perfectly elastic-plastic. Creep at contact can be enabled.\n\n.. note::\n This law uses :yref:`ScGeom`.",
((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,momentRotationLaw,false,"use bending/twisting moment at contacts. See :yref:`CohFrictPhys::cohesionDisablesFriction` for details."))
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-17 17:37:54 +0000
@@ -10,7 +10,7 @@
/********************** Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM ****************************/
CREATE_LOGGER(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM);
-void Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* scene){
+void Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
const Real& dt = scene->dt;
ScGeom* geom = static_cast<ScGeom*>(ig.get());
@@ -92,7 +92,7 @@
Vector3r f = phys->normalForce + shearForce;
// these lines to adapt to periodic boundary conditions (NOTE applyForceAtContactPoint computes torque induced by normal and shear force too)
if (!scene->isPeriodic)
- applyForceAtContactPoint(f , geom->contactPoint , id2, st2->se3.position, id1, st1->se3.position, scene);
+ applyForceAtContactPoint(f , geom->contactPoint , id2, st2->se3.position, id1, st1->se3.position);
else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used when scene is periodic
scene->forces.addForce(id1,-f);
scene->forces.addForce(id2,f);
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp 2010-07-12 17:36:42 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.hpp 2010-07-17 17:37:54 +0000
@@ -104,7 +104,7 @@
/** 2d functor creating the interaction law (Law2) based on SphereContactGeometry (ScGeom) and CFpmPhys of 2 bodies, returning type CohesiveFrictionalPM */
class Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM: public LawFunctor{
public:
- virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
FUNCTOR2D(ScGeom,CFpmPhys);
YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CFpmPhys_CohesiveFrictionalPM,LawFunctor,"Constitutive law for the CFpm model.",
=== modified file 'pkg/dem/Engine/GlobalEngine/CundallStrack.cpp'
--- pkg/dem/Engine/GlobalEngine/CundallStrack.cpp 2010-05-03 12:17:44 +0000
+++ pkg/dem/Engine/GlobalEngine/CundallStrack.cpp 2010-07-17 17:37:54 +0000
@@ -12,13 +12,13 @@
/********************** Law2_Dem3DofGeom_RockPMPhys_Rpm ****************************/
CREATE_LOGGER(Law2_Dem3Dof_CSPhys_CundallStrack);
-void Law2_Dem3Dof_CSPhys_CundallStrack::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* rootBody){
+void Law2_Dem3Dof_CSPhys_CundallStrack::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
CSPhys* phys=static_cast<CSPhys*>(ip.get());
/*NormalForce */
Real displN=geom->displacementN();
- if (displN>0){rootBody->interactions->requestErase(contact->getId1(),contact->getId2()); return;}
+ if (displN>0){scene->interactions->requestErase(contact->getId1(),contact->getId2()); return;}
phys->normalForce=phys->kn*displN*geom->normal;
/*ShearForce*/
@@ -28,7 +28,7 @@
trialFs*=sqrt(maxFsSq/(trialFs.squaredNorm()));}
phys->shearForce = trialFs;
- applyForceAtContactPoint(phys->normalForce + trialFs, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position, rootBody);
+ applyForceAtContactPoint(phys->normalForce + trialFs, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position);
return;
}
=== modified file 'pkg/dem/Engine/GlobalEngine/CundallStrack.hpp'
--- pkg/dem/Engine/GlobalEngine/CundallStrack.hpp 2010-03-15 14:35:44 +0000
+++ pkg/dem/Engine/GlobalEngine/CundallStrack.hpp 2010-07-17 17:37:54 +0000
@@ -14,7 +14,7 @@
class Law2_Dem3Dof_CSPhys_CundallStrack: public LawFunctor{
public:
- virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
FUNCTOR2D(Dem3DofGeom,CSPhys);
YADE_CLASS_BASE_DOC(Law2_Dem3Dof_CSPhys_CundallStrack,LawFunctor,"Basic constitutive law published originally by Cundall&Strack; it has normal and shear stiffnesses (Kn, Kn) and dry Coulomb friction. Operates on associated :yref:`Dem3DofGeom` and :yref:`CSPhys` instances.");
DECLARE_LOGGER;
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-17 17:37:54 +0000
@@ -30,7 +30,7 @@
// these checks would be redundant in the functor (LawDispatcher does that already)
if(!dynamic_cast<ScGeom*>(I->interactionGeometry.get()) || !dynamic_cast<FrictPhys*>(I->interactionPhysics.get())) continue;
#endif
- functor->go(I->interactionGeometry, I->interactionPhysics, I.get(), scene);
+ functor->go(I->interactionGeometry, I->interactionPhysics, I.get());
}
}
@@ -47,7 +47,7 @@
}
CREATE_LOGGER(Law2_ScGeom_FrictPhys_Basic);
-void Law2_ScGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb){
+void Law2_ScGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
const Real& dt = scene->dt;
int id1 = contact->getId1(), id2 = contact->getId2();
@@ -57,10 +57,10 @@
if (neverErase) {
currentContactPhysics->shearForce = Vector3r::Zero();
currentContactPhysics->normalForce = Vector3r::Zero();}
- else ncb->interactions->requestErase(id1,id2);
+ else scene->interactions->requestErase(id1,id2);
return;}
- State* de1 = Body::byId(id1,ncb)->state.get();
- State* de2 = Body::byId(id2,ncb)->state.get();
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
Real& un=currentContactGeometry->penetrationDepth;
TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
@@ -102,20 +102,20 @@
}
}
if (!scene->isPeriodic)
- applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position, ncb);
+ applyForceAtContactPoint(-currentContactPhysics->normalForce-shearForce, currentContactGeometry->contactPoint, id1, de1->se3.position, id2, de2->se3.position);
else {//we need to use correct branches in the periodic case, the following apply for spheres only
Vector3r force = -currentContactPhysics->normalForce-shearForce;
- ncb->forces.addForce(id1,force);
- ncb->forces.addForce(id2,-force);
- ncb->forces.addTorque(id1,(currentContactGeometry->radius1-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
- ncb->forces.addTorque(id2,(currentContactGeometry->radius2-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
+ scene->forces.addForce(id1,force);
+ scene->forces.addForce(id2,-force);
+ scene->forces.addTorque(id1,(currentContactGeometry->radius1-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
+ scene->forces.addTorque(id2,(currentContactGeometry->radius2-0.5*currentContactGeometry->penetrationDepth)* currentContactGeometry->normal.cross(force));
}
//FIXME : make sure currentContactPhysics->prevNormal is not used anywhere, remove it from physics and remove this line :
currentContactPhysics->prevNormal = currentContactGeometry->normal;
}
// same as elasticContactLaw, but using Dem3DofGeom
-void Law2_Dem3DofGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene*){
+void Law2_Dem3DofGeom_FrictPhys_Basic::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
FrictPhys* phys=static_cast<FrictPhys*>(ip.get());
Real displN=geom->displacementN();
@@ -133,6 +133,6 @@
if(trialFs.squaredNorm()>maxFsSq){ geom->slipToDisplacementTMax(sqrt(maxFsSq)/phys->ks); trialFs*=sqrt(maxFsSq/(trialFs.squaredNorm()));}
//Workaround end
phys->shearForce=trialFs;
- applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position,scene);
+ applyForceAtContactPoint(phys->normalForce+trialFs,geom->contactPoint,contact->getId1(),geom->se31.position,contact->getId2(),geom->se32.position);
}
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-07-08 18:14:12 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-07-17 17:37:54 +0000
@@ -18,7 +18,7 @@
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*);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
Real elasticEnergy ();
Real getPlasticDissipation();
void initPlasticDissipation(Real initVal=0);
@@ -44,7 +44,7 @@
*/
class Law2_Dem3DofGeom_FrictPhys_Basic: public LawFunctor{
public:
- virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene*);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
FUNCTOR2D(Dem3DofGeom,FrictPhys);
YADE_CLASS_BASE_DOC(Law2_Dem3DofGeom_FrictPhys_Basic,LawFunctor,"Constitutive law for linear compression, no tension, and linear plasticity surface.\n\nThis class serves also as tutorial and is documented in detail at https://yade-dem.org/index.php/ConstitutiveLawHowto.";);
};
=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-07-17 17:37:54 +0000
@@ -66,14 +66,14 @@
/******************** Law2_ScGeom_MindlinPhys_Mindlin *********/
CREATE_LOGGER(Law2_ScGeom_MindlinPhys_Mindlin);
-void Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* ncb){
+void Law2_ScGeom_MindlinPhys_Mindlin::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
const Real& dt = scene->dt; // get time step
int id1 = contact->getId1(); // get id body 1
int id2 = contact->getId2(); // get id body 2
- State* de1 = Body::byId(id1,ncb)->state.get();
- State* de2 = Body::byId(id2,ncb)->state.get();
+ State* de1 = Body::byId(id1,scene)->state.get();
+ State* de2 = Body::byId(id2,scene)->state.get();
ScGeom* scg = static_cast<ScGeom*>(ig.get());
MindlinPhys* phys = static_cast<MindlinPhys*>(ip.get());
@@ -81,7 +81,7 @@
/*** NORMAL FORCE ***/
Real uN = scg->penetrationDepth; // get overlapping
- if (uN<0) {ncb->interactions->requestErase(contact->getId1(),contact->getId2()); return;}
+ if (uN<0) {scene->interactions->requestErase(contact->getId1(),contact->getId2()); return;}
/* 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)
In the first place we get the normal force and then we store kn to be passed to the GSTS */
@@ -140,13 +140,13 @@
/*** APPLY FORCES ***/
if (!scene->isPeriodic)
- applyForceAtContactPoint(-phys->normalForce - phys->shearForce, scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position, ncb);
+ applyForceAtContactPoint(-phys->normalForce - phys->shearForce, scg->contactPoint , id1, de1->se3.position, id2, de2->se3.position);
else { // in scg we do not wrap particles positions, hence "applyForceAtContactPoint" cannot be used
Vector3r force = -phys->normalForce - phys->shearForce;
- ncb->forces.addForce(id1,force);
- ncb->forces.addForce(id2,-force);
- ncb->forces.addTorque(id1,(scg->radius1-0.5*scg->penetrationDepth)* scg->normal.cross(force));
- ncb->forces.addTorque(id2,(scg->radius2-0.5*scg->penetrationDepth)* scg->normal.cross(force));
+ scene->forces.addForce(id1,force);
+ scene->forces.addForce(id2,-force);
+ scene->forces.addTorque(id1,(scg->radius1-0.5*scg->penetrationDepth)* scg->normal.cross(force));
+ scene->forces.addTorque(id2,(scg->radius2-0.5*scg->penetrationDepth)* scg->normal.cross(force));
}
phys->prevNormal = scg->normal;
}
=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp 2010-06-21 17:31:31 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.hpp 2010-07-17 17:37:54 +0000
@@ -51,7 +51,7 @@
class Law2_ScGeom_MindlinPhys_Mindlin: public LawFunctor{
public:
Real cn, cs;
- virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
FUNCTOR2D(ScGeom,MindlinPhys);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_MindlinPhys_Mindlin,LawFunctor,"Constitutive law for the Hertz-Mindlin formulation. It includes non linear elasticity in the normal direction as predicted by Hertz for two non-conforming elastic contact bodies. In the shear direction, instead, it reseambles the simplified case without slip discussed in Mindlin's paper, where a linear relationship between shear force and tangential displacement is provided. Finally, the Mohr-Coulomb criterion is employed to established the maximum friction force which can be developed at the contact. Moreover, it is also possible to include the effect of contact damping through the definition of the parameters $\\beta_{n}$ and $\\beta_{s}$.",
((bool,preventGranularRatcheting,false,"bool to avoid granular ratcheting"))
=== modified file 'pkg/dem/meta/ConcretePM.cpp'
--- pkg/dem/meta/ConcretePM.cpp 2010-06-29 14:32:03 +0000
+++ pkg/dem/meta/ConcretePM.cpp 2010-07-17 17:37:54 +0000
@@ -4,7 +4,7 @@
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/pkg-dem/Shop.hpp>
-YADE_PLUGIN((CpmState)(CpmMat)(Ip2_CpmMat_CpmMat_CpmPhys)(CpmPhys)(Law2_Dem3DofGeom_CpmPhys_Cpm)
+YADE_PLUGIN((CpmState)(CpmMat)(Ip2_CpmMat_CpmMat_CpmPhys)(CpmPhys)(Law2_Dem3DofGeom_CpmPhys_Cpm)(Law2_ScGeom_CpmPhys_Cpm)
#ifdef YADE_OPENGL
(Gl1_CpmPhys)
#endif
@@ -137,7 +137,7 @@
#endif
-void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* scene){
+void Law2_Dem3DofGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I){
Dem3DofGeom* contGeom=static_cast<Dem3DofGeom*>(_geom.get());
CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
@@ -158,22 +158,19 @@
Real& epsNPl(BC->epsNPl); const Real& dt=scene->dt; const Real& dmgTau(BC->dmgTau); const Real& plTau(BC->plTau);const Real& yieldLogSpeed(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldLogSpeed); const int& yieldSurfType(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSurfType); const Real& yieldEllipseShift(Law2_Dem3DofGeom_CpmPhys_Cpm::yieldEllipseShift); const Real& epsSoft(Law2_Dem3DofGeom_CpmPhys_Cpm::epsSoft); const Real& relKnSoft(Law2_Dem3DofGeom_CpmPhys_Cpm::relKnSoft);
#endif
-
epsN=contGeom->strainN(); epsT=contGeom->strainT();
// debugging
- #define YADE_VERIFY(condition) if(!(condition)){LOG_FATAL("Verification `"<<#condition<<"' failed!"); LOG_FATAL("in interaction #"<<I->getId1()<<"+#"<<I->getId2()); Omega::instance().saveSimulation("/tmp/verificationFailed.xml"); throw;}
-
- #define NNAN(a) YADE_VERIFY(!isnan(a));
- #define NNANV(v) YADE_VERIFY(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
- #ifdef YADE_DEBUG
- if(isnan(epsN)){
- LOG_FATAL("refLength="<<contGeom->refLength<<"; pos1="<<contGeom->se31.position<<"; pos2="<<contGeom->se32.position<<"; displacementN="<<contGeom->displacementN());
- throw runtime_error("!! epsN==NaN !!");
- }
- NNAN(epsN); NNANV(epsT);
- #endif
- NNAN(epsN); NNANV(epsT);
+ #define YADE_VERIFY(condition) if(!(condition)){LOG_FATAL("Verification `"<<#condition<<"' failed!"); LOG_FATAL("in interaction #"<<I->getId1()<<"+#"<<I->getId2()); Omega::instance().saveSimulation("/tmp/verificationFailed.xml"); throw;}
+ #define NNAN(a) YADE_VERIFY(!isnan(a));
+ #define NNANV(v) YADE_VERIFY(!isnan(v[0])); assert(!isnan(v[1])); assert(!isnan(v[2]));
+ #ifdef YADE_DEBUG
+ if(isnan(epsN)){
+ LOG_FATAL("refLength="<<contGeom->refLength<<"; pos1="<<contGeom->se31.position<<"; pos2="<<contGeom->se32.position<<"; displacementN="<<contGeom->displacementN());
+ throw runtime_error("!! epsN==NaN !!");
+ }
+ #endif
+ NNAN(epsN); NNANV(epsT);
// constitutive law
#ifdef CPM_MATERIAL_MODEL
@@ -186,7 +183,6 @@
else epsN+=sigmaSoft/E+(BC->isoPrestress-sigmaSoft)/(E*relKnSoft);
}
CPM_MATERIAL_MODEL
- //cerr<<"kappaD="<<kappaD<<", epsN="<<epsN<<", neverDamage="<<neverDamage<<", epsCrackOnset="<<epsCrackOnset<<", epsFracture="<<epsFracture<<", omega="<<omega<<"="<<funcG(kappaD,epsCrackOnset,epsFracture,neverDamage);
#else
// simplified public model
epsN+=BC->isoPrestress/E;
@@ -223,7 +219,8 @@
Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
Fs=sigmaT*crossSection; BC->shearForce=Fs;
- applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(), contGeom->se32.position, scene);
+ applyForceAtContactPoint(BC->normalForce+BC->shearForce, contGeom->contactPoint, I->getId1(), contGeom->se31.position, I->getId2(), contGeom->se32.position);
+ #undef YADE_VERIFY
}
Real Law2_Dem3DofGeom_CpmPhys_Cpm::yieldSigmaTMagnitude(Real sigmaN, Real omega, Real undamagedCohesion, Real tanFrictionAngle){
@@ -236,6 +233,69 @@
}
+CREATE_LOGGER(Law2_ScGeom_CpmPhys_Cpm);
+void Law2_ScGeom_CpmPhys_Cpm::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I){
+ ScGeom* contGeom=static_cast<ScGeom*>(_geom.get());
+ CpmPhys* BC=static_cast<CpmPhys*>(_phys.get());
+ // FIXME: better way to get current positions? Needed for
+ // (1) getting intial equilibrium distance (not stored in ScGeom and distFactor not accessible from here)
+ // (2) applying contact forces back on particles
+ const Vector3r& pos1(Body::byId(I->getId1(),scene)->state->pos); const Vector3r& pos2(Body::byId(I->getId2(),scene)->state->pos);
+ // just the first time
+ if(I->isFresh(scene)){
+ Real minRad=(contGeom->refR1<=0?contGeom->refR2:(contGeom->refR2<=0?contGeom->refR1:min(contGeom->refR1,contGeom->refR2)));
+ BC->crossSection=Mathr::PI*pow(minRad,2);
+ BC->refLength=(pos2-pos1).norm();
+ BC->kn=BC->crossSection*BC->E/BC->refLength;
+ BC->ks=BC->crossSection*BC->G/BC->refLength;
+ }
+ // shorthands
+ Real& epsN(BC->epsN);
+ Vector3r& epsT(BC->epsT); Real& kappaD(BC->kappaD); Real& epsPlSum(BC->epsPlSum); const Real& E(BC->E); const Real& undamagedCohesion(BC->undamagedCohesion); const Real& tanFrictionAngle(BC->tanFrictionAngle); const Real& G(BC->G); const Real& crossSection(BC->crossSection); const Real& omegaThreshold(Law2_ScGeom_CpmPhys_Cpm::omegaThreshold); const Real& epsCrackOnset(BC->epsCrackOnset); Real& relResidualStrength(BC->relResidualStrength); const Real& epsFracture(BC->epsFracture); const bool& neverDamage(BC->neverDamage); Real& omega(BC->omega); Real& sigmaN(BC->sigmaN); Vector3r& sigmaT(BC->sigmaT); Real& Fn(BC->Fn); Vector3r& Fs(BC->Fs); // for python access
+ const bool& isCohesive(BC->isCohesive);
+ epsN=contGeom->penetrationDepth/BC->refLength;
+
+ // FIXME: sign?
+ epsT=contGeom->rotate(epsT);
+ epsT+=contGeom->shearIncrement()/BC->refLength;
+
+ // simplified public model
+ epsN+=BC->isoPrestress/E;
+ // very simplified version of the constitutive law
+ kappaD=max(max(0.,epsN),kappaD); // internal variable, max positive strain (non-decreasing)
+ omega=isCohesive?Law2_Dem3DofGeom_CpmPhys_Cpm::funcG(kappaD,epsCrackOnset,epsFracture,neverDamage):1.; // damage variable (non-decreasing, as funcG is also non-decreasing)
+ sigmaN=(1-(epsN>0?omega:0))*E*epsN; // damage taken in account in tension only
+ sigmaT=G*epsT; // trial stress
+ Real yieldSigmaT=max((Real)0.,undamagedCohesion*(1-omega)-sigmaN*tanFrictionAngle); // Mohr-Coulomb law with damage
+ if(sigmaT.squaredNorm()>yieldSigmaT*yieldSigmaT){
+ Real scale=yieldSigmaT/sigmaT.norm();
+ sigmaT*=scale;
+ epsPlSum+=(epsT-epsT*scale).norm()*yieldSigmaT;
+ epsT*=scale;
+ }
+ relResidualStrength=isCohesive?(kappaD<epsCrackOnset?1.:(1-omega)*(kappaD)/epsCrackOnset):0;
+ sigmaN-=BC->isoPrestress;
+
+ // handle broken contacts
+ if(epsN>0. && ((isCohesive && omega>omegaThreshold) || !isCohesive)){
+ if(isCohesive){
+ const shared_ptr<Body>& body1=Body::byId(I->getId1(),scene), body2=Body::byId(I->getId2(),scene); assert(body1); assert(body2);
+ const shared_ptr<CpmState>& st1=YADE_PTR_CAST<CpmState>(body1->state), st2=YADE_PTR_CAST<CpmState>(body2->state);
+ // nice article about openMP::critical vs. scoped locks: http://www.thinkingparallel.com/2006/08/21/scoped-locking-vs-critical-in-openmp-a-personal-shootout/
+ { 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->getId1(),I->getId2());
+ return;
+ }
+
+ Fn=sigmaN*crossSection; BC->normalForce=Fn*contGeom->normal;
+ Fs=sigmaT*crossSection; BC->shearForce=Fs;
+
+ applyForceAtContactPoint(BC->normalForce+BC->shearForce,contGeom->contactPoint,I->getId1(),pos1,I->getId2(),pos2);
+}
+
+
#ifdef YADE_OPENGL
/********************** Gl1_CpmPhys ****************************/
#include<yade/lib-opengl/OpenGLWrapper.hpp>
=== modified file 'pkg/dem/meta/ConcretePM.hpp'
--- pkg/dem/meta/ConcretePM.hpp 2010-07-12 17:36:42 +0000
+++ pkg/dem/meta/ConcretePM.hpp 2010-07-17 17:37:54 +0000
@@ -109,7 +109,7 @@
static long cummBetaIter, cummBetaCount;
/*! auxiliary variable for visualization, recalculated in Law2_Dem3DofGeom_CpmPhys_Cpm at every iteration */
// Fn and Fs are also stored as Vector3r normalForce, shearForce in NormShearPhys
- Real omega, Fn, sigmaN, epsN, relResidualStrength; Vector3r epsT, sigmaT, Fs;
+ Real omega, Fn, sigmaN, epsN, relResidualStrength; Vector3r sigmaT, Fs;
static Real solveBeta(const Real c, const Real N);
Real computeDmgOverstress(Real dt);
@@ -137,11 +137,12 @@
((Real,epsTrans,0,"Transversal strain (perpendicular to the contact axis)"))
((Real,epsPlSum,0,"cummulative shear plastic strain measure (scalar) on this contact"))
((bool,isCohesive,false,"if not cohesive, interaction is deleted when distance is greater than zero."))
+ ((Vector3r,epsT,Vector3r::Zero(),"Total shear strain (either computed from increments with :yref:`ScGeom` or simple copied with :yref:`Dem3DofGeom`) |yupdate|"))
+ ((Real,refLength,NaN,"Reference contact length (only used with :yref:`Law2_ScGeom_CpmPhys_Cpm`)"))
,
createIndex(); epsT=Fs=Vector3r::Zero(); Fn=0; omega=0;
,
.def_readonly("omega",&CpmPhys::omega,"Damage internal variable")
- .def_readonly("epsT",&CpmPhys::epsT,"Transversal strain (not used)")
.def_readonly("Fn",&CpmPhys::Fn,"Magnitude of normal force.")
.def_readonly("Fs",&CpmPhys::Fs,"Magnitude of shear force")
.def_readonly("epsN",&CpmPhys::epsN,"Current normal strain")
@@ -176,14 +177,14 @@
class Law2_Dem3DofGeom_CpmPhys_Cpm: public LawFunctor{
public:
/*! Damage evolution law */
- Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) {
+ static Real funcG(const Real& kappaD, const Real& epsCrackOnset, const Real& epsFracture, const bool& neverDamage) {
if(kappaD<epsCrackOnset || neverDamage) return 0;
return 1.-(epsCrackOnset/kappaD)*exp(-(kappaD-epsCrackOnset)/epsFracture);
}
//! return |sigmaT| at plastic surface for given sigmaN etc; not used by the law itself
Real yieldSigmaTMagnitude(Real sigmaN, Real omega, Real undamagedCohesion, Real tanFrictionAngle);
- void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+ void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
FUNCTOR2D(Dem3DofGeom,CpmPhys);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_Dem3DofGeom_CpmPhys_Cpm,LawFunctor,"Constitutive law for the :ref:`cpm-model`.",
@@ -201,6 +202,20 @@
};
REGISTER_SERIALIZABLE(Law2_Dem3DofGeom_CpmPhys_Cpm);
+class Law2_ScGeom_CpmPhys_Cpm: public LawFunctor{
+ public:
+ void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
+ FUNCTOR2D(ScGeom,CpmPhys);
+ YADE_CLASS_BASE_DOC_ATTRS(Law2_ScGeom_CpmPhys_Cpm,LawFunctor,"An experimental version of :yref:`Law2_Dem3DofGeom_CpmPhys_Cpm` which uses :yref:`ScGeom` instead of :yref:`Dem3DofGeom`.",
+ ((Real,omegaThreshold,((void)">=1. to deactivate, i.e. never delete any contacts",1.),"damage after which the contact disappears (<1), since omega reaches 1 only for strain →+∞"))
+ );
+ DECLARE_LOGGER;
+};
+REGISTER_SERIALIZABLE(Law2_ScGeom_CpmPhys_Cpm);
+
+
+
+
#ifdef YADE_OPENGL
#include<yade/pkg-common/GLDrawFunctors.hpp>
class Gl1_CpmPhys: public GlInteractionPhysicsFunctor {
=== modified file 'pkg/dem/meta/RockPM.cpp'
--- pkg/dem/meta/RockPM.cpp 2010-07-14 11:30:24 +0000
+++ pkg/dem/meta/RockPM.cpp 2010-07-17 17:37:54 +0000
@@ -10,7 +10,7 @@
/********************** Law2_Dem3DofGeom_RockPMPhys_Rpm ****************************/
CREATE_LOGGER(Law2_Dem3DofGeom_RockPMPhys_Rpm);
-void Law2_Dem3DofGeom_RockPMPhys_Rpm::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact, Scene* rootBody){
+void Law2_Dem3DofGeom_RockPMPhys_Rpm::go(shared_ptr<InteractionGeometry>& ig, shared_ptr<InteractionPhysics>& ip, Interaction* contact){
Dem3DofGeom* geom=static_cast<Dem3DofGeom*>(ig.get());
RpmPhys* phys=static_cast<RpmPhys*>(ip.get());
@@ -52,7 +52,7 @@
phys->shearForce = Fs;
/**Normal Interaction*/
- applyForceAtContactPoint(phys->normalForce + phys->shearForce, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position, rootBody);
+ applyForceAtContactPoint(phys->normalForce + phys->shearForce, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position);
if ((phys->isCohesive)&&(displN<(-phys->lengthMaxCompression))) {
phys->isCohesive = false; ///Destruction
@@ -67,7 +67,7 @@
* Destruction.
**/
if (displN>(phys->lengthMaxTension)) {
- rootBody->interactions->requestErase(contact->getId1(),contact->getId2());
+ scene->interactions->requestErase(contact->getId1(),contact->getId2());
return;
} else {
/**If the distance
@@ -75,14 +75,14 @@
* we aply additional forces to keep particles together.
**/
phys->normalForce=phys->kn*displN*geom->normal;
- applyForceAtContactPoint(phys->normalForce, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position, rootBody);
+ applyForceAtContactPoint(phys->normalForce, geom->contactPoint, contact->getId1(), geom->se31.position, contact->getId2(), geom->se32.position);
return;
}
} else {
/**
* Delete interactions
*/
- rootBody->interactions->requestErase(contact->getId1(),contact->getId2());
+ scene->interactions->requestErase(contact->getId1(),contact->getId2());
return;
}
}
=== modified file 'pkg/dem/meta/RockPM.hpp'
--- pkg/dem/meta/RockPM.hpp 2010-07-14 12:44:25 +0000
+++ pkg/dem/meta/RockPM.hpp 2010-07-17 17:37:54 +0000
@@ -28,7 +28,7 @@
class Law2_Dem3DofGeom_RockPMPhys_Rpm: public LawFunctor{
public:
- virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody);
+ virtual void go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I);
FUNCTOR2D(Dem3DofGeom,RpmPhys);
YADE_CLASS_BASE_DOC(Law2_Dem3DofGeom_RockPMPhys_Rpm,LawFunctor,"Constitutive law for the Rpm model");
=== modified file 'pkg/dem/meta/Shop.cpp'
--- pkg/dem/meta/Shop.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/meta/Shop.cpp 2010-07-17 17:37:54 +0000
@@ -135,11 +135,11 @@
/* Apply force on contact point to 2 bodies; the force is oriented as it applies on the first body and is reversed on the second.
*/
-void Shop::applyForceAtContactPoint(const Vector3r& force, const Vector3r& contPt, body_id_t id1, const Vector3r& pos1, body_id_t id2, const Vector3r& pos2, Scene* rootBody){
- rootBody->forces.addForce(id1,force);
- rootBody->forces.addForce(id2,-force);
- rootBody->forces.addTorque(id1,(contPt-pos1).cross(force));
- rootBody->forces.addTorque(id2,-(contPt-pos2).cross(force));
+void Shop::applyForceAtContactPoint(const Vector3r& force, const Vector3r& contPt, body_id_t id1, const Vector3r& pos1, body_id_t id2, const Vector3r& pos2, Scene* scene){
+ scene->forces.addForce(id1,force);
+ scene->forces.addForce(id2,-force);
+ scene->forces.addTorque(id1,(contPt-pos1).cross(force));
+ scene->forces.addTorque(id2,-(contPt-pos2).cross(force));
}
=== modified file 'pkg/dem/meta/Shop.hpp'
--- pkg/dem/meta/Shop.hpp 2010-05-19 11:09:23 +0000
+++ pkg/dem/meta/Shop.hpp 2010-07-17 17:37:54 +0000
@@ -108,7 +108,7 @@
static shared_ptr<Interaction> createExplicitInteraction(body_id_t id1, body_id_t id2, bool force);
//! apply force on contact point on both bodies (reversed on body 2)
- static void applyForceAtContactPoint(const Vector3r& force, const Vector3r& contPt, body_id_t id1, const Vector3r& pos1, body_id_t id2, const Vector3r& pos2, Scene* rb);
+ static void applyForceAtContactPoint(const Vector3r& force, const Vector3r& contPt, body_id_t id1, const Vector3r& pos1, body_id_t id2, const Vector3r& pos2, Scene* scene);
//! map scalar variable to 1d colorscale
static Vector3r scalarOnColorScale(Real x, Real xmin=0., Real xmax=1.);
=== modified file 'pkg/dem/meta/ViscoelasticPM.cpp'
--- pkg/dem/meta/ViscoelasticPM.cpp 2010-05-13 20:19:38 +0000
+++ pkg/dem/meta/ViscoelasticPM.cpp 2010-07-17 17:37:54 +0000
@@ -32,7 +32,7 @@
}
/* Law2_ScGeom_ViscElPhys_Basic */
-void Law2_ScGeom_ViscElPhys_Basic::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I, Scene* rootBody){
+void Law2_ScGeom_ViscElPhys_Basic::go(shared_ptr<InteractionGeometry>& _geom, shared_ptr<InteractionPhysics>& _phys, Interaction* I){
const ScGeom& geom=*static_cast<ScGeom*>(_geom.get());
ViscElPhys& phys=*static_cast<ViscElPhys*>(_phys.get());
@@ -41,18 +41,18 @@
const int id2 = I->getId2();
if (geom.penetrationDepth<0) {
- rootBody->interactions->requestErase(id1,id2);
+ scene->interactions->requestErase(id1,id2);
return;
}
- const BodyContainer& bodies = *rootBody->bodies;
+ const BodyContainer& bodies = *scene->bodies;
const State& de1 = *static_cast<State*>(bodies[id1]->state.get());
const State& de2 = *static_cast<State*>(bodies[id2]->state.get());
Vector3r& shearForce = phys.shearForce;
- if (I->isFresh(rootBody)) shearForce=Vector3r(0,0,0);
+ if (I->isFresh(scene)) shearForce=Vector3r(0,0,0);
const Real& dt = scene->dt;
@@ -99,9 +99,9 @@
}
const Vector3r f = phys.normalForce + shearForce + shearForceVisc;
- addForce (id1,-f,rootBody);
- addForce (id2, f,rootBody);
- addTorque(id1,-c1x.cross(f),rootBody);
- addTorque(id2, c2x.cross(f),rootBody);
+ addForce (id1,-f,scene);
+ addForce (id2, f,scene);
+ addTorque(id1,-c1x.cross(f),scene);
+ addTorque(id2, c2x.cross(f),scene);
}
=== modified file 'pkg/dem/meta/ViscoelasticPM.hpp'
--- pkg/dem/meta/ViscoelasticPM.hpp 2010-04-10 15:11:48 +0000
+++ pkg/dem/meta/ViscoelasticPM.hpp 2010-07-17 17:37:54 +0000
@@ -57,7 +57,7 @@
/// This class provides linear viscoelastic contact model
class Law2_ScGeom_ViscElPhys_Basic: public LawFunctor {
public :
- virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*, Scene*);
+ virtual void go(shared_ptr<InteractionGeometry>&, shared_ptr<InteractionPhysics>&, Interaction*);
FUNCTOR2D(ScGeom,ViscElPhys);
YADE_CLASS_BASE_DOC(Law2_ScGeom_ViscElPhys_Basic,LawFunctor,"Linear viscoelastic model operating on :yref:`ScGeom` and :yref:`ViscElPhys`.");
};
=== modified file 'py/plot.py'
--- py/plot.py 2010-07-08 17:19:27 +0000
+++ py/plot.py 2010-07-17 17:37:54 +0000
@@ -164,6 +164,8 @@
pylab.xlabel(xlateLabel(p))
if 'title' in O.tags.keys(): pylab.title(O.tags['title'])
+
+
def liveUpdate(timestamp):
global liveTimeStamp
liveTimeStamp=timestamp
@@ -208,7 +210,18 @@
thread.start_new_thread(liveUpdate,(time.time(),))
# pylab.show() # this blocks for some reason; call show on figures directly
figs=set([l.line.get_axes().get_figure() for l in currLineRefs])
- for f in figs: f.show()
+ for f in figs:
+ f.show()
+ # should have fixed https://bugs.launchpad.net/yade/+bug/606220, but does not work apparently
+ if 0:
+ import matplotlib.backend_bases
+ if 'CloseEvent' in dir(matplotlib.backend_bases):
+ def closeFigureCallback(event):
+ ff=event.canvas.figure
+ # remove closed axes from our update list
+ global currLineRefs
+ currLineRefs=[l for l in currLineRefs if l.line.get_axes().get_figure()!=ff]
+ f.canvas.mpl_connect('close_event',closeFigureCallback)
else: return pylab.gcf() # return the current figure
def saveGnuplot(baseName,term='wxt',extension=None,timestamp=False,comment=None,title=None,varData=False):
Follow ups
