yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05434
[Branch ~yade-dev/yade/trunk] Rev 2377: Restore the part of r2367 concerning ScGeom (+ some more improvments).
------------------------------------------------------------
revno: 2377
committer: bchareyre <bchareyre@dt-rv020>
branch nick: trunk
timestamp: Sun 2010-07-18 21:21:08 +0200
message:
Restore the part of r2367 concerning ScGeom (+ some more improvments).
modified:
pkg/common/DataClass/Shape/Cylinder.cpp
pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp
pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp
pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp
pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp
pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp
pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp
py/utils.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/DataClass/Shape/Cylinder.cpp'
--- pkg/common/DataClass/Shape/Cylinder.cpp 2010-07-15 14:26:48 +0000
+++ pkg/common/DataClass/Shape/Cylinder.cpp 2010-07-18 19:21:08 +0000
@@ -52,18 +52,13 @@
Real norm=normal.norm(); normal/=norm; // normal is unit vector now
Real penetrationDepth=s1->radius+s2->radius-norm;
scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
- if(isNew) scm->prevNormal=normal;
- else scm->prevNormal=scm->normal;
- scm->normal=normal;
+// if(isNew) scm->prevNormal=normal;
+// else scm->prevNormal=scm->normal;
+// scm->normal=normal;
scm->penetrationDepth=penetrationDepth;
scm->radius1=s1->radius;
scm->radius2=s2->radius;
-#ifdef IGCACHE
-// if (scene->isPeriodic) {
-// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
-// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
- /*else */scm->precompute(state1,state2,scene,c,true);
-#endif
+ scm->precompute(state1,state2,scene,c,normal,isNew,true);
return true;
}
return false;
@@ -114,26 +109,20 @@
else { scm=shared_ptr<ScGeom>(new ScGeom()); c->interactionGeometry=scm; }
Real length=(bchain2.pos-bchain1.pos).norm();
Vector3r segt =pChain2->pos-pChain1->pos;
- if(isNew) {scm->normal=scm->prevNormal=segt/length;bs1->initLength=length;}
- else {scm->prevNormal=scm->normal; scm->normal=segt/length;}
+ if(isNew) {/*scm->normal=scm->prevNormal=segt/length;*/bs1->initLength=length;}
scm->radius1=revert ? 0:bs1->initLength;
scm->radius2=revert ? bs1->initLength:0;
scm->penetrationDepth=bs1->initLength-length;
scm->contactPoint=bchain2.pos;
//bs1->segment used for fast BBs and projections + display
- if (revert) segt = -segt;
- bs1->segment=segt;
+ bs1->segment= bchain2.pos-bchain1.pos;
#ifdef YADE_OPENGL
//bs1->length and s1->chainedOrientation used for display only,
bs1->length=length;
bs1->chainedOrientation.setFromTwoVectors(Vector3r::UnitZ(),bchain1.ori.conjugate()*segt);
#endif
+ scm->precompute(state1,state2,scene,c,segt/length,isNew,true);
-#ifdef IGCACHE
- scm->precompute(state1,state2,scene,c,true);
-#else
- cerr<<"this is not supposed to work without #define IGCACHE"<<endl;
-#endif
//Set values that will be considered in Ip2 functor, geometry (precomputed) is really defined with values above
scm->radius1 = scm->radius2 = bs1->initLength*0.5;
return true;
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-07-18 19:21:08 +0000
@@ -7,10 +7,7 @@
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-
YADE_PLUGIN((ScGeom));
-
-// At least one virtual method must be in the .cpp file (!!!)
ScGeom::~ScGeom(){};
Vector3r& ScGeom::rotate(Vector3r& shearForce) const {
@@ -23,13 +20,22 @@
}
//!Precompute data needed for rotating tangent vectors attached to the interaction
-void ScGeom::precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, bool avoidGranularRatcheting){
- orthonormal_axis = prevNormal.cross(normal);
- Real angle = scene->dt*0.5*normal.dot(rbp1.angVel + rbp2.angVel);
- twist_axis = angle*normal;
- //The previous normal could be updated here after being used, them remove it from Ig2's?
- //prevNormal=normal;
+void ScGeom::precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r& currentNormal, bool isNew, bool avoidGranularRatcheting){
+ if(!isNew) {
+ orthonormal_axis = normal.cross(currentNormal);
+ Real angle = scene->dt*0.5*normal.dot(rbp1.angVel + rbp2.angVel);
+ twist_axis = angle*normal;}
+ else twist_axis=orthonormal_axis=Vector3r::Zero();
+ //Update contact normal
+ normal=currentNormal;
//Precompute shear increment
+ Vector3r relativeVelocity=getIncidentVel(&rbp1,&rbp2,scene->dt,scene->isPeriodic ? Vector3r(scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>()) : Vector3r::Zero(),avoidGranularRatcheting);
+ //keep the shear part only
+ relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
+ shearInc = relativeVelocity*scene->dt;
+}
+
+Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
Vector3r& x = contactPoint;
Vector3r c1x, c2x;
if(avoidGranularRatcheting){
@@ -54,93 +60,6 @@
c2x = -radius2*normal;
} else {
// FIXME: It is correct for sphere-sphere and sphere-facet contact
- c1x = (x - rbp1.pos);
- c2x = (x - rbp2.pos);
- }
- Vector3r relativeVelocity = (rbp2.vel+rbp2.angVel.cross(c2x)) - (rbp1.vel+rbp1.angVel.cross(c1x));
- if (scene->isPeriodic) relativeVelocity += scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
- //keep the shear part only
- relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
- shearInc = relativeVelocity*scene->dt;
-}
-
-
-//Do we need a function for such simple operation? BTW, the original was BAD duplicated code. ;-)
-Vector3r ScGeom::updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){rotate(shear); shear -= shearInc;}
-
-#ifndef IGCACHE
-Vector3r ScGeom::rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
- //Just pass null shift to the periodic version
- return rotateAndGetShear(shearForce,prevNormal,rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
-}
-#endif
-
-//Removing this function will need to adapt all laws
-const Vector3r& ScGeom::rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
-#ifdef IGCACHE
- rotate(shearForce);
- return shearInc;
- #else
- rotate(shearForce,prevNormal,rbp1,rbp2,dt);
- Vector3r& x = contactPoint;
- Vector3r c1x, c2x;
- if(avoidGranularRatcheting){
- // FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
- c1x = radius1*normal;
- c2x = -radius2*normal;
- } else {
- // FIXME: It is correct for sphere-sphere and sphere-facet contact
- c1x = (x - rbp1->pos);
- c2x = (x - rbp2->pos);
- }
- Vector3r relativeVelocity = (rbp2->vel+rbp2->angVel.cross(c2x)) - (rbp1->vel+rbp1->angVel.cross(c1x));
-
- //update relative velocity for interactions across periods
- relativeVelocity+=shiftVel;
- //keep the shear part only
- relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
- Vector3r shearDisplacement = relativeVelocity*dt;
- return shearDisplacement;
-#endif //IGCACHE
-}
-
-Vector3r& ScGeom::rotate(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt){
-#ifdef IGCACHE
- return rotate(shearForce);
-#else
- Vector3r axis;
- Real angle;
- // approximated rotations
- axis = prevNormal.cross(normal);
- shearForce -= shearForce.cross(axis);
- angle = dt*0.5*normal.dot(rbp1->angVel + rbp2->angVel);
- axis = angle*normal;
- shearForce -= shearForce.cross(axis);
- // exact rotations
-
-// Quaternionr q;
-// axis = prevNormal.Cross(normal);
-// angle = acos(normal.Dot(prevNormal));
-// q.FromAngleAxis(angle,axis);
-// shearForce = shearForce*q;
-// angle = dt*0.5*normal.dot(rbp1->angVel+rbp2->angVel);
-// axis = normal;
-// q.FromAngleAxis(angle,axis);
-// shearForce = q*shearForce;
-
- return shearForce;
-#endif
-}
-
-Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
- Vector3r& x = contactPoint;
- Vector3r c1x, c2x;
- if(avoidGranularRatcheting){
- // FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
- c1x = radius1*normal;
- c2x = -radius2*normal;
- } else {
- // FIXME: It is correct for sphere-sphere and sphere-facet contact
c1x = (x - rbp1->pos);
c2x = (x - rbp2->pos);
}
@@ -154,22 +73,3 @@
//Just pass null shift to the periodic version
return getIncidentVel(rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
}
-
-/* keep this for reference; declarations in the header */
-#if 0
- Vector3r ScGeom::relRotVector() const{
- Quaternionr relOri12=ori1.Conjugate()*ori2;
- Quaternionr oriDiff=initRelOri12.Conjugate()*relOri12;
- Vector3r axis; Real angle;
- oriDiff.ToAxisAngle(axis,angle);
- if(angle>Mathr::PI)angle-=Mathr::TWO_PI;
- return angle*axis;
- }
-
- void ScGeom::bendingTorsionAbs(Vector3r& bend, Real& tors){
- Vector3r relRot=relRotVector();
- tors=relRot.Dot(normal);
- bend=relRot-tors*normal;
- }
-#endif
-
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-07-18 19:21:08 +0000
@@ -28,23 +28,16 @@
Real penetrationDepth;
virtual ~ScGeom();
- //!precompute values of shear increment and interaction rotation data
- void precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, bool avoidGranularRatcheting=true);
- //! Periodic variant. Needs the velocity shift across periods for periodic BCs (else it is safe to pass Vector3r::Zero()). Typically obtained as scene->cell->velGrad*scene->cell->Hsize*cellDist.
-// void precompute(const State& rbp1, const State& rbp2, const Real& dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
+ //!precompute values of shear increment and interaction rotation data. Update contact normal to the vurrentNormal value. Precondition : the value of normal is not updated outside (and before) this function.
+ void precompute(const State& rbp1, const State& rbp2, const Scene* scene, const shared_ptr<Interaction>& c, const Vector3r& currentNormal, bool isNew, bool avoidGranularRatcheting=true);
//! Rotates a "shear" vector to keep track of contact orientation. Returns reference of the updated vector.
Vector3r& rotate(Vector3r& tangentVector) const;
const Vector3r& shearIncrement() const {return shearInc;}
//!________________REMOVED IN NEXT STEP___________________________________
- Vector3r& rotate(Vector3r& tangentVector, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt);
//! update shear on this contact given dynamic parameters of bodies. Should be called from constitutive law, exactly once per iteration. Returns shear increment in this update, which is already added to shear.
- Vector3r updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
- //! Update the shear force orientation (or any other vector, think "shear") on this contact given dynamic parameters of bodies. Compute and return shear displacement increment.
- const Vector3r& rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt,bool avoidGranularRatcheting=true);
- //! Periodic variant. Needs the velocity shift across periods for periodic BCs (else it is safe to pass Vector3r::Zero()). Typically obtained as scene->cell->velGrad*scene->cell->Hsize*cellDist. It would be better to define the shift transparently inside the function, but it needs scene and interaction pointers, which we don't have here.
- const Vector3r& rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
+ void updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
//!________________END REMOVED IN NEXT STEP___________________________________
// Add method which returns the impact velocity (then, inside the contact law, this can be split into shear and normal component). Handle periodicity.
@@ -53,9 +46,7 @@
Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting=true);
YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`bodies<Body>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
- ((Vector3r,contactPoint,Vector3r::Zero(),"Reference point of the contact. |ycomp|"))
- ((Vector3r,shear,Vector3r::Zero(),"Total value of the current shear. Update the value using ScGeom::updateShear. This is experimental. You have to update it yourself in constitutive laws using e.g. :yref:`ScGeom::rotate` and :yref:`ScGeom::shearIncrement`"))
- ((Vector3r,prevNormal,Vector3r::Zero(),"Normal of the contact in the previous step. |ycomp|")) ,
+ ((Vector3r,contactPoint,Vector3r::Zero(),"Reference point of the contact. |ycomp|")),
/* extra initializers */ ((radius1,GenericSpheresContact::refR1)) ((radius2,GenericSpheresContact::refR2)),
/* ctor */ createIndex(); shearInc=Vector3r::Zero(); twist_axis=orthonormal_axis=Vector3r::Zero();,
/* py */
=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-07-18 19:21:08 +0000
@@ -99,18 +99,16 @@
bool isNew=!c->interactionGeometry;
if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
-
- if(isNew) { /* same as below */ scm->prevNormal=normal;}
- else {scm->prevNormal=scm->normal;}
// contact point is in the middle of overlapping volumes
//(in the direction of penetration, which is normal to the box surface closest to sphere center) of overlapping volumes
scm->contactPoint = 0.5*(pt1+pt2);
- scm->normal = normal;
+// scm->normal = normal;
scm->penetrationDepth = (pt1-pt2).norm();
scm->radius1 = s->radius;
scm->radius2 = s->radius;
c->interactionGeometry = scm;
+ scm->precompute(state1,state2,scene,c,normal,isNew,true);
} else { // outside
Vector3r cOnBox_box = boxAxisT*cOnBox_boxLocal; // projection of sphere's center on closest box surface - relative to box's origin, but GLOBAL orientation!
Vector3r cOnBox_sphere = cOnBox_box-relPos21; // same, but origin in sphere's center
@@ -144,23 +142,13 @@
bool isNew=!c->interactionGeometry;
if (isNew) scm = shared_ptr<ScGeom>(new ScGeom());
else scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
- if(isNew) { /* same as below */ scm->prevNormal=-cOnBox_sphere; }
- else scm->prevNormal=scm->normal;
scm->contactPoint = 0.5*(pt1+pt2);
- //scm->normal = pt1-pt2; scm->normal.normalize();
- //scm->penetrationDepth = (pt1-pt2).norm();
- scm->normal = -cOnBox_sphere;
scm->penetrationDepth = depth;
scm->radius1 = s->radius;
scm->radius2 = s->radius;
c->interactionGeometry = scm;
+ scm->precompute(state1,state2,scene,c,-cOnBox_sphere,isNew,true);
}
-#ifdef IGCACHE
-// if (scene->isPeriodic){
-// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
-// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
- /*else */scm->precompute(state1,state2,scene,c,true);
-#endif
return true;
}
=== modified file 'pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp 2010-07-18 19:21:08 +0000
@@ -103,19 +103,11 @@
normal = facetAxisT*normal; // in global orientation
scm->contactPoint = se32.position - (sphereRadius-0.5*penetrationDepth)*normal;
- //Update prevNormal (mandatory in ScGeom algoritms)
- if(isNew) { scm->prevNormal=normal;}
- else scm->prevNormal=scm->normal;
- scm->normal = normal;
scm->penetrationDepth = penetrationDepth;
scm->radius1 = 2*sphereRadius;
scm->radius2 = sphereRadius;
- if (!c->interactionGeometry)
- c->interactionGeometry = scm;
-// if (scene->isPeriodic){
-// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
-// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
- /*else */scm->precompute(state1,state2,scene,c,true);
+ if (isNew) c->interactionGeometry = scm;
+ scm->precompute(state1,state2,scene,c,normal,isNew,true);
return true;
}
return false;
=== modified file 'pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-07-15 14:26:48 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-07-18 19:21:08 +0000
@@ -36,23 +36,10 @@
Real norm=normal.norm(); normal/=norm; // normal is unit vector now
Real penetrationDepth=s1->radius+s2->radius-norm;
scm->contactPoint=se31.position+(s1->radius-0.5*penetrationDepth)*normal;//0.5*(pt1+pt2);
- if(isNew) scm->prevNormal=normal;
- else scm->prevNormal=scm->normal;
- scm->normal=normal;
scm->penetrationDepth=penetrationDepth;
scm->radius1=s1->radius;
scm->radius2=s2->radius;
-#ifdef IGCACHE
-// if (scene->isPeriodic) {
-// Vector3r shiftVel = scene->cell->velGrad*scene->cell->Hsize*c->cellDist.cast<Real>();
-// scm->precompute(state1,state2,scene->dt,shiftVel,true);}
- /*else */scm->precompute(state1,state2,scene,c,true);
-#endif
-
- // keep this for reference on how to compute bending and torsion from relative orientation; parts in ScGeom header
- #if 0
- scm->initRelOri12=se31.orientation.Conjugate()*se32.orientation;
- #endif
+ scm->precompute(state1,state2,scene,c,normal,isNew,true);
return true;
}
return false;
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-17 17:37:54 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-18 19:21:08 +0000
@@ -74,12 +74,10 @@
///////////////////////// CREEP START ///////////
if (shear_creep) shearForce -= currentContactPhysics->ks*(shearForce*dt/creep_viscosity);
///////////////////////// CREEP END ////////////
-#ifdef IGCACHE
+
Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
const Vector3r& dus = currentContactGeometry->shearIncrement();
-#else
- Vector3r dus = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt);
-#endif
+
//Linear elasticity giving "trial" shear force
shearForce -= currentContactPhysics->ks*dus;
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-18 11:59:57 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-18 19:21:08 +0000
@@ -65,15 +65,9 @@
State* st1 = Body::byId(id1,scene)->state.get();
State* st2 = Body::byId(id2,scene)->state.get();
-#ifdef IGCACHE
geom->rotate(phys->shearForce);
const Vector3r& dus = geom->shearIncrement();
-#else
- // define shift to handle periodicity
- Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
- Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, shiftVel, preventGranularRatcheting);
-#endif
- /// before changes to adapt periodic scene it was like that: Vector3r dus = geom->rotateAndGetShear(shearForce, phys->prevNormal, st1, st2, dt, preventGranularRatcheting);
+
//Linear elasticity giving "trial" shear force
shearForce -= phys->ks*dus;
// needed for the next timestep
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-17 17:37:54 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-18 19:21:08 +0000
@@ -21,7 +21,6 @@
void ElasticContactLaw::action()
{
if(!functor) functor=shared_ptr<Law2_ScGeom_FrictPhys_Basic>(new Law2_ScGeom_FrictPhys_Basic);
- functor->useShear=useShear;
functor->neverErase=neverErase;
functor->scene=scene;
FOREACH(const shared_ptr<Interaction>& I, *scene->interactions){
@@ -48,7 +47,6 @@
CREATE_LOGGER(Law2_ScGeom_FrictPhys_Basic);
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();
ScGeom* currentContactGeometry= static_cast<ScGeom*>(ig.get());
@@ -64,30 +62,20 @@
Real& un=currentContactGeometry->penetrationDepth;
TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
-#ifdef IGCACHE
+
Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
const Vector3r& shearDisp = currentContactGeometry->shearIncrement();
-#else
- Vector3r& shearForce = currentContactPhysics->shearForce;
- Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
- Vector3r shearDisp = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt,shiftVel,true);
-// cerr << "shearForce "<<shearForce<<" shearIncrement "<<currentContactGeometry->shearIncrement<<endl;
-#endif
+
if (!traceEnergy){//Update force but don't compute energy terms (see below))
- if(useShear){
- currentContactGeometry->updateShear(de1,de2,dt);
- shearForce=currentContactPhysics->ks*currentContactGeometry->shear;
- } else shearForce -= currentContactPhysics->ks*shearDisp;
+ shearForce -= currentContactPhysics->ks*shearDisp;
// PFC3d SlipModel, is using friction angle. CoulombCriterion
Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
if( shearForce.squaredNorm() > maxFs ){
Real ratio = Mathr::Sqrt(maxFs) / shearForce.norm();
- shearForce *= ratio;
- if(useShear) currentContactGeometry->shear*=ratio;}
+ shearForce *= ratio;}
} else {
- //almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false
- if(useShear) throw ("energy tracing not defined with useShear==true");
+ //almost the same with additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vector when traceEnergy==false
shearForce -= currentContactPhysics->ks*shearDisp;
Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-07-17 17:37:54 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-07-18 19:21:08 +0000
@@ -24,7 +24,6 @@
void initPlasticDissipation(Real initVal=0);
YADE_CLASS_BASE_DOC_ATTRS_CTOR_PY(Law2_ScGeom_FrictPhys_Basic,LawFunctor,"Law for linear compression, without cohesion and Mohr-Coulomb plasticity surface.\n\n.. note::\n This law uses :yref:`ScGeom`; there is also functionally equivalent :yref:`Law2_Dem3DofGeom_FrictPhys_Basic`, which uses :yref:`Dem3DofGeom` (sphere-box interactions are not implemented for the latest).",
((bool,neverErase,false,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
- ((bool,useShear,false,"Use ScGeom::updateShear rather than ScGeom::rotateAndGetShear for shear force computation."))
((bool,traceEnergy,false,"Define the total energy dissipated in plastic slips at all contacts."))
,,
.def("elasticEnergy",&Law2_ScGeom_FrictPhys_Basic::elasticEnergy,"Compute and return the total elastic energy in all \"FrictPhys\" contacts")
@@ -56,7 +55,6 @@
void action();
YADE_CLASS_BASE_DOC_ATTRS(ElasticContactLaw,GlobalEngine,"[DEPRECATED] Loop over interactions applying :yref:`Law2_ScGeom_FrictPhys_Basic` on all interactions.\n\n.. note::\n Use :yref:`InteractionDispatchers` and :yref:`Law2_ScGeom_FrictPhys_Basic` instead of this class for performance reasons.",
((bool,neverErase,false,"Keep interactions even if particles go away from each other (only in case another constitutive law is in the scene, e.g. :yref:`Law2_ScGeom_CapillaryPhys_Capillarity`)"))
- ((bool,useShear,false,"Use :yref:`ScGeom`::updateShear rather than :yref:`ScGeom`::rotateAndGetShear for shear force computation."))
);
};
=== modified file 'pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp'
--- pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-07-17 17:37:54 +0000
+++ pkg/dem/Engine/GlobalEngine/HertzMindlin.cpp 2010-07-18 19:21:08 +0000
@@ -112,7 +112,8 @@
// Define shift to handle periodicity
Vector3r shiftVel = scene->isPeriodic ? (Vector3r)((scene->cell->velGrad*scene->cell->Hsize)*Vector3r((Real) contact->cellDist[0],(Real) contact->cellDist[1],(Real) contact->cellDist[2])) : Vector3r::Zero();
// 1. Rotate shear force
- shearElastic = scg->rotate(shearElastic, phys->prevNormal, de1, de2, dt);
+// shearElastic = scg->rotate(shearElastic, phys->prevNormal, de1, de2, dt);
+ shearElastic = scg->rotate(shearElastic);
// 2. Get incident velocity, get shear and normal components
// FIXME: Concerning the possibility to avoid granular ratcheting, it is not clear how this works in the case of HM. See the thread http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg01947.html
Vector3r incidentV = scg->getIncidentVel(de1, de2, dt, shiftVel, preventGranularRatcheting);
=== modified file 'py/utils.py'
--- py/utils.py 2010-07-15 14:26:48 +0000
+++ py/utils.py 2010-07-18 19:21:08 +0000
@@ -238,7 +238,7 @@
segment=end-begin
b=Body()
b.shape=ChainedCylinder(radius=radius,length=segment.norm(),color=color if color else randomColor(),wire=wire,highlight=highlight)
- b.shape.segment=b.shape.length*Vector3(0.,0.,1.)
+ b.shape.segment=segment
V=2*(4./3)*math.pi*radius**3
geomInert=(2./5.)*V*radius**2+b.shape.length*b.shape.length*2*(4./3)*math.pi*radius**3
b.state=ChainedState()
Follow ups
