yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05334
[Branch ~yade-dev/yade/trunk] Rev 2360: - Cylinder : fix segment orientation, use correct contact point, and clean code.
------------------------------------------------------------
revno: 2360
committer: bchareyre <bchareyre@dt-rv020>
branch nick: trunk
timestamp: Tue 2010-07-13 21:14:55 +0200
message:
- Cylinder : fix segment orientation, use correct contact point, and clean code.
- ScGeom : last step before removing old code, some functors adapted (others are using even older duplicated code apparently)
- Small sphinx fixes.
modified:
doc/references.bib
doc/sphinx/formulation.rst
doc/yade-articles.bib
pkg/common/DataClass/Shape/Cylinder.cpp
pkg/common/DataClass/Shape/Cylinder.hpp
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
scripts/test/chained-cylinder-spring.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 'doc/references.bib'
--- doc/references.bib 2010-07-12 14:26:27 +0000
+++ doc/references.bib 2010-07-13 19:14:55 +0000
@@ -73,18 +73,31 @@
volume = {131},
number = {7},
pages = {689-698},
- url = {https://yade-dem.org/wiki/File:Chareyre%26Villard2005_licensed.pdf},
+ url = {http://link.aip.org/link/?QEM/131/689/1},
doi = {10.1061/(ASCE)0733-9399(2005)131:7(689)}
}
@phdthesis{Scholtes2009,
- author = {Luc Scholtès},
- title = {Modélisation micromécanique des milieux granulaires partiellement saturés},
+ author = { Luc Scholtès },
+ title = { Modélisation micromécanique des milieux granulaires partiellement saturés },
year = {2009},
school = {Institut National Polytechnique de Grenoble},
url = {http://tel.archives-ouvertes.fr/tel-00363961/en/}
}
+@article{Scholtes2009b,
+ author = "L. Scholtès and B. Chareyre and F. Nicot and F. Darve",
+ title = "Micromechanics of granular materials with capillary effects",
+ journal = "International Journal of Engineering Science",
+ volume = "47",
+ number = "1",
+ pages = "64--75",
+ year = "2009",
+ issn = "0020-7225",
+ doi = "10.1016/j.ijengsci.2008.07.002",
+ url = "http://www.sciencedirect.com/science/article/B6V32-4TDJGDJ-1/2/9878efa5e573197aff6ee14d5abb58a2";,
+}
+
@proceedings{DeghmReport2006,
title={Annual Report 2006},
year={2006},
@@ -94,6 +107,18 @@
url={ http://geo.hmg.inpg.fr/frederic/Discrete_Element_Group_FVD.html }
}
+@article{Jerier2009,
+ journal={Granular Matter},
+ title={A geometric algorithm based on tetrahedral meshes to generate a dense polydisperse sphere packing},
+ volume={11},
+ number={1},
+ pages={43--52},
+ year={2009},
+ author={Jean-François Jerier and Didier Imbault and Fréderic-Victor Donzé and Pierre Doremus},
+ doi={10.1007/s10035-008-0116-0},
+ url={http://www.springerlink.com/content/w0x307g110421035}
+}
+
@article{Pournin2001,
title = {Molecular-dynamics force models for better control of energy dissipation in numerical simulations of dense granular media},
author = {Pournin, L. and Liebling, Th. M. and Mocellin, A. },
=== modified file 'doc/sphinx/formulation.rst'
--- doc/sphinx/formulation.rst 2010-07-08 20:01:59 +0000
+++ doc/sphinx/formulation.rst 2010-07-13 19:14:55 +0000
@@ -717,7 +717,7 @@
DEM simulations
""""""""""""""""
-In DEM simulations, per-particle stiffness $\vec{K}_{ij}$ is determined from the stiffnesses of contacts in which it participates [Chareyre2005]_. Suppose each contact has normal stiffness $K_{Nk}$, shear stiffness $K_{Tk}=\xi K_{Nk}$ and is oriented by normal $\vec{n}_{k}$. A translational stiffness matrix $\vec{K}_ij$ can be defined as the sum of contributions of all contacts in which it participates (indices $k$), as
+In DEM simulations, per-particle stiffness $\vec{K}_{ij}$ is determined from the stiffnesses of contacts in which it participates [Chareyre2005]_. Suppose each contact has normal stiffness $K_{Nk}$, shear stiffness $K_{Tk}=\xi K_{Nk}$ and is oriented by normal $\vec{n}_{k}$. A translational stiffness matrix $\vec{K}_{ij}$ can be defined as the sum of contributions of all contacts in which it participates (indices $k$), as
.. math::
:label: eq-dtcr-particle-stiffness
=== modified file 'doc/yade-articles.bib'
--- doc/yade-articles.bib 2010-07-05 08:16:58 +0000
+++ doc/yade-articles.bib 2010-07-13 19:14:55 +0000
@@ -149,7 +149,8 @@
journal = { J. Engrg. Mech.},
volume = {131},
number = {7},
- year = {2005}
+ year = {2005},
+ url={https://yade-dem.org/wiki/File:Chareyre%26Villard2005_licensed.pdf}
}
@article{Hentz2004a,
=== modified file 'pkg/common/DataClass/Shape/Cylinder.cpp'
--- pkg/common/DataClass/Shape/Cylinder.cpp 2010-07-12 19:41:57 +0000
+++ pkg/common/DataClass/Shape/Cylinder.cpp 2010-07-13 19:14:55 +0000
@@ -16,7 +16,7 @@
#ifdef YADE_OPENGL
(Gl1_Cylinder)(Gl1_ChainedCylinder)
#endif
- (Bo1_Cylinder_Aabb)/*(Bo1_ChainedCylinder_Aabb)*/
+ (Bo1_Cylinder_Aabb)(Bo1_ChainedCylinder_Aabb)
);
vector<vector<int> > ChainedState::chains;
@@ -59,10 +59,10 @@
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->dt,true);
+// 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;
}
@@ -97,65 +97,48 @@
const ChainedState *pChain1, *pChain2;
pChain1=YADE_CAST<const ChainedState*>(&state1);
pChain2=YADE_CAST<const ChainedState*>(&state2);
- const ChainedState& chain1 = *pChain1;
- const ChainedState& chain2 = *pChain2;
if (!pChain1 || !pChain2) {
cerr <<"cast failed8567"<<endl;
}
- bool revert = (pChain2->rank-chain1.rank == -1);
-
+ const bool revert = (pChain2->rank-pChain1->rank == -1);
+ const ChainedState& bchain1 = revert? *pChain2 : *YADE_CAST<const ChainedState*>(&state1);
+ const ChainedState& bchain2 = revert? *pChain1 : *pChain2;
+ if (bchain2.rank-bchain1.rank != 1) {/*cerr<<"Mutual contacts in same chain between not adjacent elements, not handled*/ return false;}
if (pChain2->chainNumber!=pChain1->chainNumber) {cerr<<"PROBLEM0124"<<endl; return false;}
- else if (chain2.rank-chain1.rank != 1 && pChain2->rank-pChain1->rank != -1) {/*cerr<<"PROBLEM0123"<<endl;*/ return false;}
- const ChainedState* pChain = revert ? pChain2 : pChain1;
-
- ChainedCylinder *s1=static_cast<ChainedCylinder*>(cm1.get()), *s2=static_cast<ChainedCylinder*>(cm2.get());
- ChainedCylinder *s = revert ? s2 : s1;
+
+ ChainedCylinder *bs1=static_cast<ChainedCylinder*>(revert? cm2.get():cm1.get());
+
shared_ptr<ScGeom> scm;
bool isNew = !c->interactionGeometry;
if(!isNew) scm=YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
else { scm=shared_ptr<ScGeom>(new ScGeom()); c->interactionGeometry=scm; }
-// Real length=abs((pChain->ori*Vector3r::UnitZ()).dot(chain2.pos-chain1.pos));
- Real length=(chain2.pos-chain1.pos).norm();
-// Vector3r segment =pChain->ori*Vector3r::UnitZ()*length;
- Vector3r segment =chain2.pos-chain1.pos;
- if (revert) segment = -segment;
- s->chainedOrientation.setFromTwoVectors(Vector3r::UnitZ(),pChain->ori.conjugate()*segment);
-// s->chainedOrientation=pChain->ori.conjugate()*s->chainedOrientation;
-
- if (revert) segment = -segment;
- if(isNew) {scm->normal=scm->prevNormal=segment/length;s->initLength=length;}
- else {scm->prevNormal=scm->normal; scm->normal=segment/length;}
-
- if (!revert) {
- scm->radius1=s->initLength;
- scm->radius2=0;}
- else {
- scm->radius1=0;
- scm->radius2=s->initLength;}
- //length only used for display
- s->length=length;
- scm->penetrationDepth=s->initLength-length;
-
-// scm->contactPoint=pChain->pos+segment;//generates instabilites (TODO)
- scm->contactPoint=pChain->pos+pChain->ori*Vector3r::UnitZ()*length;
+ 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;}
+ 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;
+#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
#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->dt,true);
- }
+ 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=s->initLength*0.5;
- scm->radius2=s->initLength*0.5;
+ scm->radius1 = scm->radius2 = bs1->initLength*0.5;
return true;
}
-
bool Ig2_ChainedCylinder_ChainedCylinder_ScGeom::goReverse( const shared_ptr<Shape>& cm1,
const shared_ptr<Shape>& cm2,
const State& state1,
@@ -176,12 +159,6 @@
int Gl1_Cylinder::glutStacks;
int Gl1_Cylinder::glCylinderList=-1;
-// bool Gl1_ChainedCylinder::wire;
-// bool Gl1_ChainedCylinder::glutNormalize;
-// int Gl1_ChainedCylinder::glutSlices;
-// int Gl1_ChainedCylinder::glutStacks;
-// int Gl1_ChainedCylinder::glCylinderList=-1;
-
void Gl1_Cylinder::out( Quaternionr q )
{
AngleAxisr aa(q);
@@ -194,7 +171,7 @@
Real length=(static_cast<Cylinder*>(cm.get()))->length;
//glMaterialv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, Vector3f(cm->color[0],cm->color[1],cm->color[2]));
glColor3v(cm->color);
- if(glutNormalize) glPushAttrib(GL_NORMALIZE); // as per http://lists.apple.com/archives/Mac-opengl/2002/Jul/msg00085.html
+ if(glutNormalize) glPushAttrib(GL_NORMALIZE);
// glPushMatrix();
Quaternionr shift = (static_cast<ChainedCylinder*>(cm.get()))->chainedOrientation;
if (wire || wire2) drawCylinder(true, r,length,shift);
@@ -204,20 +181,17 @@
return;
}
-void Gl1_ChainedCylinder::go(const shared_ptr<Shape>& cm, const shared_ptr<State>& ,bool wire2, const GLViewInfo&)
+void Gl1_ChainedCylinder::go(const shared_ptr<Shape>& cm, const shared_ptr<State>& state,bool wire2, const GLViewInfo&)
{
Real r=(static_cast<ChainedCylinder*>(cm.get()))->radius;
Real length=(static_cast<ChainedCylinder*>(cm.get()))->length;
- Quaternionr shift = (static_cast<ChainedCylinder*>(cm.get()))->chainedOrientation;
- //glMaterialv(GL_FRONT, GL_AMBIENT_AND_DIFFUSE, Vector3f(cm->color[0],cm->color[1],cm->color[2]));
+ Quaternionr shift;// = (static_cast<ChainedCylinder*>(cm.get()))->chainedOrientation;
+ shift.setFromTwoVectors(Vector3r::UnitZ(),state->ori.conjugate()*(static_cast<ChainedCylinder*>(cm.get()))->segment);
glColor3v(cm->color);
- if(glutNormalize) glPushAttrib(GL_NORMALIZE); // as per http://lists.apple.com/archives/Mac-opengl/2002/Jul/msg00085.html
-// glPushMatrix();
- //out(shift);
+ if(glutNormalize) glPushAttrib(GL_NORMALIZE);
if (wire || wire2) drawCylinder(true, r,length,shift);
else drawCylinder(false, r,length,shift);
if(glutNormalize) glPopAttrib();
-// glPopMatrix();
return;
}
@@ -316,15 +290,14 @@
// aabb->min = scene->cell->unshearPt(se3.position)-minkSize;
// aabb->max = scene->cell->unshearPt(se3.position)+minkSize;
}
-/*
+
+
void Bo1_ChainedCylinder_Aabb::go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r& se3, const Body* b){
Cylinder* cylinder = static_cast<Cylinder*>(cm.get());
Aabb* aabb = static_cast<Aabb*>(bv.get());
if(!scene->isPeriodic){
- const Vector3r& O = se3.position-se3.positioncylinder->segment*0.5;
- //In this case, segment will be updated by Ig2_ChainedCylinder_ChainedCylinder_ScGeom, no need to rotate it
- Vector3r O2 = se3.position+cylinder->segment*0.5;
-// cerr << O<<" "<<O2<<endl;
+ const Vector3r& O = se3.position;
+ Vector3r O2 = se3.position+cylinder->segment;
aabb->min=aabb->max=O;
for (int k=0;k<3;k++){
aabb->min[k]=min(aabb->min[k],min(O[k],O2[k])-cylinder->radius);
@@ -332,8 +305,5 @@
}
return;
}
-}*/
-
-
-
+}
=== modified file 'pkg/common/DataClass/Shape/Cylinder.hpp'
--- pkg/common/DataClass/Shape/Cylinder.hpp 2010-07-12 19:41:57 +0000
+++ pkg/common/DataClass/Shape/Cylinder.hpp 2010-07-13 19:14:55 +0000
@@ -4,6 +4,7 @@
#include<yade/core/Body.hpp>
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/InteractionGeometryFunctor.hpp>
+#include<yade/pkg-common/Sphere.hpp>
#include<yade/core/Scene.hpp>
#ifdef YADE_OPENGL
#include<yade/pkg-common/GLDrawFunctors.hpp>
@@ -11,17 +12,18 @@
#include<yade/pkg-common/BoundFunctor.hpp>
-class Cylinder: public Shape{
+class Cylinder: public Sphere{
public:
- Cylinder(Real _radius, Real _length): radius(_radius), length(_length), segment(Vector3r(0,0,1)*_length){ createIndex(); }
+ Cylinder(Real _radius, Real _length): length(_length) { /*segment=Vector3r(0,0,1)*_length;*/ radius=_radius; createIndex(); }
virtual ~Cylinder ();
- YADE_CLASS_BASE_DOC_ATTRS_CTOR(Cylinder,Shape,"Geometry of a cylinder, as Minkowski sum of line and sphere.",
- ((Real,radius,NaN,"Radius [m]"))
+ YADE_CLASS_BASE_DOC_ATTRS_CTOR(Cylinder,Sphere,"Geometry of a cylinder, as Minkowski sum of line and sphere.",
+// ((Real,radius,NaN,"Radius [m]"))
((Real,length,NaN,"Length [m]"))
((Vector3r,segment,Vector3r::Zero(),"Length vector")),
- createIndex(); /*ctor*/
+ createIndex();
+ /*ctor*/
);
- REGISTER_CLASS_INDEX(Cylinder,Shape);
+ REGISTER_CLASS_INDEX(Cylinder,Sphere);
};
class ChainedCylinder: public Cylinder{
@@ -142,7 +144,7 @@
//!This doesn't work : the 1D dispatcher will pick Gl1_Cylinder to display ChainedCylinders, workaround : add shift to cylinders (should be a variable of chained cylinders only).
class Gl1_ChainedCylinder : public Gl1_Cylinder{
public:
- virtual void go(const shared_ptr<Shape>&, const shared_ptr<State>&,bool,const GLViewInfo&);
+ virtual void go(const shared_ptr<Shape>&, const shared_ptr<State>& state, bool,const GLViewInfo&);
YADE_CLASS_BASE_DOC(Gl1_ChainedCylinder,Gl1_Cylinder,"Renders :yref:`ChainedCylinder` object including a shift for compensating flexion."
);
RENDERS(ChainedCylinder);
@@ -180,6 +182,18 @@
);
};
+class Bo1_ChainedCylinder_Aabb : public BoundFunctor
+{
+ public :
+ void go(const shared_ptr<Shape>& cm, shared_ptr<Bound>& bv, const Se3r&, const Body*);
+ FUNCTOR2D(ChainedCylinder,Aabb);
+ YADE_CLASS_BASE_DOC_ATTRS(Bo1_ChainedCylinder_Aabb,BoundFunctor,"Functor creating :yref:`Aabb` from :yref:`ChainedCylinder`.",
+ ((Real,aabbEnlargeFactor,((void)"deactivated",-1),"Relative enlargement of the bounding box; deactivated if negative.\n\n.. note::\n\tThis attribute is used to create distant interaction, but is only meaningful with an :yref:`InteractionGeometryFunctor` which will not simply discard such interactions: :yref:`Ig2_Cylinder_Cylinder_Dem3DofGeom::distFactor` / :yref:`Ig2_Cylinder_Cylinder_ScGeom::interactionDetectionFactor` should have the same value as :yref:`aabbEnlargeFactor<Bo1_Cylinder_Aabb::aabbEnlargeFactor>`."))
+ );
+};
+
+
+
// Keep this : Cylinders and ChainedCylinders will have different centers maybe.
// class Bo1_ChainedCylinder_Aabb : public Bo1_Cylinder_Aabb
// {
@@ -195,7 +209,7 @@
REGISTER_SERIALIZABLE(Bo1_Cylinder_Aabb);
-// REGISTER_SERIALIZABLE(Bo1_ChainedCylinder_Aabb);
+REGISTER_SERIALIZABLE(Bo1_ChainedCylinder_Aabb);
#ifdef YADE_OPENGL
REGISTER_SERIALIZABLE(Gl1_Cylinder);
REGISTER_SERIALIZABLE(Gl1_ChainedCylinder);
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-07-12 08:50:35 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-07-13 19:14:55 +0000
@@ -5,12 +5,15 @@
#include "ScGeom.hpp"
#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){
+Vector3r& ScGeom::rotate(Vector3r& shearForce) const {
// approximated rotations
shearForce -= shearForce.cross(orthonormal_axis);
shearForce -= shearForce.cross(twist_axis);
@@ -19,15 +22,10 @@
return shearForce;
}
-void ScGeom::precompute(const State& rbp1, const State& rbp2, const Real& dt,bool avoidGranularRatcheting){
- //Pass null shift to the periodic version
- precompute(rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
-}
-
-void ScGeom::precompute(const State& rbp1, const State& rbp2, const Real& dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
-//Precompute data needed for rotating tangent vectors attached to the interaction
+//!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 = dt*0.5*normal.dot(rbp1.angVel + rbp2.angVel);
+ 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;
@@ -60,76 +58,29 @@
c2x = (x - rbp2.pos);
}
Vector3r relativeVelocity = (rbp2.vel+rbp2.angVel.cross(c2x)) - (rbp1.vel+rbp1.angVel.cross(c1x));
- relativeVelocity+=shiftVel;
+ 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;
- shearIncrement = relativeVelocity*dt;
-}
-
-
-
-Vector3r ScGeom::updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
-
- Vector3r axis;
- Real angle;
-
- Vector3r shearIncrement(Vector3r::Zero());
-
- // approximated rotations
- axis = prevNormal.cross(normal);
- shearIncrement -= shear.cross(axis);
- angle = dt*0.5*normal.dot(rbp1->angVel + rbp2->angVel);
- axis = angle*normal;
- shearIncrement -= (shear+shearIncrement).cross(axis);
-
- // exact rotations (not adapted to shear/shearIncrement!)
- #if 0
- 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;
- #endif
-
- 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));
- Vector3r shearVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
- Vector3r shearDisplacement = shearVelocity*dt;
- shearIncrement -= shearDisplacement;
-
- shear+=shearIncrement;
- return shearIncrement;
-}
-
-//Removing this function will need to adapt all laws
+ 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
-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 shearIncrement;
-// #else
+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;
@@ -150,23 +101,23 @@
relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
Vector3r shearDisplacement = relativeVelocity*dt;
return shearDisplacement;
-// #endif //IGCACHE
+#endif //IGCACHE
}
Vector3r& ScGeom::rotate(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt){
-// #ifdef IGCACHE
-// return rotate(shearForce);
-// #else
+#ifdef IGCACHE
+ return rotate(shearForce);
+#else
Vector3r axis;
Real angle;
// approximated rotations
- axis = prevNormal.cross(normal);
+ 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));
@@ -178,7 +129,7 @@
// shearForce = q*shearForce;
return shearForce;
-// #endif
+#endif
}
Vector3r ScGeom::getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting){
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-07-12 08:50:35 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-07-13 19:14:55 +0000
@@ -3,12 +3,12 @@
// © 2006 Bruno Chareyre <bruno.chareyre@xxxxxxxxxxx>
#pragma once
-
+#include<yade/core/Interaction.hpp>
#include<yade/core/InteractionGeometry.hpp>
#include<yade/core/State.hpp>
#include<yade/lib-base/Math.hpp>
#include<yade/pkg-dem/DemXDofGeom.hpp>
-/*! Class representing geometry of two spheres in contact.
+/*! Class representing geometry of two bodies in contact.
*
* The code under SCG_SHEAR is experimental and is used only if ElasticContactLaw::useShear is explicitly true
*/
@@ -17,48 +17,50 @@
#define IGCACHE
class ScGeom: public GenericSpheresContact {
+ private:
+ //cached values
+ Vector3r twist_axis;//rotation vector arounf normal
+ Vector3r orthonormal_axis;//rotation vector in contact plane
+ Vector3r shearInc;//shear disp. increment
public:
// inherited from GenericSpheresContact: Vector3r& normal;
Real &radius1, &radius2;
- //cached values
- Vector3r twist_axis;//rotation vector arounf normal
- Vector3r orthonormal_axis;//rotation vector in contact plane
Real penetrationDepth;
- Vector3r shearIncrement;
+ 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);
+
+ //! 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);
-
- virtual ~ScGeom();
//! 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.
- Vector3r rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt,bool avoidGranularRatcheting=true);
+ 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.
- Vector3r rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
-
- //!precompute values of shear increment and interaction rotation data
- void precompute(const State& rbp1, const State& rbp2, const 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.
- void precompute(const State& rbp1, const State& rbp2, const Real& dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
-
- //! Add method which only rotates the shear vector (or another tangent vector as well). Returns reference of the updated vector.
- Vector3r& rotate(Vector3r& tangentVector, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt);
- //! Same rotation using precomputed values
- Vector3r& rotate(Vector3r& tangentVector);
+ const Vector3r& rotateAndGetShear(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, 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.
Vector3r getIncidentVel(const State* rbp1, const State* rbp2, Real dt, const Vector3r& shiftVel, bool avoidGranularRatcheting=true);
// Implement another version of getIncidentVel which does not handle periodicity.
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:`spheres<Sphere>` in contact. The contact has 3 DOFs (normal and 2×shear) and uses incremental algorithm for updating shear. (For shear formulated in total displacements and rotations, see :yref:`Dem3DofGeom` and related classes).\n\nWe use symbols $\\vec{x}$, $\\vec{v}$, $\\vec{\\omega}$ respectively for position, linear and angular velocities (all in global coordinates) and $r$ for particles radii; subscripted with 1 or 2 to distinguish 2 spheres in contact. Then we compute unit contact normal\n\n.. math::\n\n\t\\vec{n}=\\frac{\\vec{x}_2-\\vec{x}_1}{||\\vec{x}_2-\\vec{x}_1||}\n\nRelative velocity of spheres is then\n\n.. math::\n\n\t\\vec{v}_{12}=(\\vec{v}_2+\\vec{\\omega}_2\\times(-r_2\\vec{n}))-(\\vec{v}_1+\\vec{\\omega}_1\\times(r_1\\vec{n}))\n\nand its shear component\n\n.. math::\n\n\t\\Delta\\vec{v}_{12}^s=\\vec{v}_{12}-(\\vec{n}\\cdot\\vec{v}_{12})\\vec{n}.\n\nTangential displacement increment over last step then reads\n\n.. math::\n\n\t\\vec{x}_{12}^s=\\Delta t \\vec{v}_{12}^s.",
+ YADE_CLASS_BASE_DOC_ATTRS_INIT_CTOR_PY(ScGeom,GenericSpheresContact,"Class representing :yref:`geometry<InteractionGeometry>` of two :yref:`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. |ycomp|"))
- ((Vector3r,prevNormal,Vector3r::Zero(),"Normal of the contact in the previous step. |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|")) ,
/* extra initializers */ ((radius1,GenericSpheresContact::refR1)) ((radius2,GenericSpheresContact::refR2)),
- /* ctor */ createIndex(); shearIncrement=Vector3r::Zero(); twist_axis=orthonormal_axis=Vector3r::Zero();,
+ /* ctor */ createIndex(); shearInc=Vector3r::Zero(); twist_axis=orthonormal_axis=Vector3r::Zero();,
/* py */
.def_readwrite("penetrationDepth",&ScGeom::penetrationDepth,"documentation")
+ .def_readonly("shearIncrement",&ScGeom::shearInc,"documentation")
);
REGISTER_CLASS_INDEX(ScGeom,GenericSpheresContact);
};
=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-07-08 18:37:01 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-07-13 19:14:55 +0000
@@ -156,10 +156,10 @@
c->interactionGeometry = scm;
}
#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->dt,true);
+// 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-08 18:14:12 +0000
+++ pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp 2010-07-13 19:14:55 +0000
@@ -112,10 +112,10 @@
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->dt,true);
+// 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);
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-08 18:14:12 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-07-13 19:14:55 +0000
@@ -43,10 +43,10 @@
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->dt,true);
+// 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
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-12 17:44:48 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalContactLaw.cpp 2010-07-13 19:14:55 +0000
@@ -63,7 +63,6 @@
if (un < 0 && (currentContactPhysics->normalForce.squaredNorm() > pow(currentContactPhysics->normalAdhesion,2)
|| currentContactPhysics->normalAdhesion==0)) {
// BREAK due to tension
- cerr <<"requesterase";
ncb->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;
@@ -77,7 +76,7 @@
///////////////////////// CREEP END ////////////
#ifdef IGCACHE
Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
- const Vector3r& dus = currentContactGeometry->shearIncrement;
+ const Vector3r& dus = currentContactGeometry->shearIncrement();
#else
Vector3r dus = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt);
#endif
=== modified file 'pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp'
--- pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-06-08 22:25:00 +0000
+++ pkg/dem/Engine/GlobalEngine/CohesiveFrictionalPM.cpp 2010-07-13 19:14:55 +0000
@@ -66,9 +66,15 @@
// using scGeom function rotateAndGetShear
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;
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-08 18:14:12 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-13 19:14:55 +0000
@@ -66,7 +66,7 @@
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
#ifdef IGCACHE
Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
- Vector3r& shearDisp = currentContactGeometry->shearIncrement;
+ 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();
=== modified file 'scripts/test/chained-cylinder-spring.py'
--- scripts/test/chained-cylinder-spring.py 2010-07-12 19:41:57 +0000
+++ scripts/test/chained-cylinder-spring.py 2010-07-13 19:14:55 +0000
@@ -6,21 +6,21 @@
from yade import utils
young=1.0e5
-poisson=10
+poisson=5
density=2.60e3
frictionAngle=radians(30)
O.materials.append(CohFrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,label='mat'))
-O.dt=9e-5
+O.dt=1e-4
O.initializers=[
## Create bounding boxes. They are needed to zoom the 3d view properly before we start the simulation.
- BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Cylinder_Aabb()])
+ BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_ChainedCylinder_Aabb()])
]
O.engines=[
ForceResetter(),
BoundDispatcher([
- Bo1_Cylinder_Aabb(),
+ Bo1_ChainedCylinder_Aabb(),
Bo1_Sphere_Aabb()
]),
InsertionSortCollider(),
@@ -34,7 +34,7 @@
## Motion equation
NewtonIntegrator(damping=0.15),
PeriodicPythonRunner(iterPeriod=500,command='history()'),
- #PeriodicPythonRunner(iterPeriod=100,command='yade.qt.center()')
+ PeriodicPythonRunner(iterPeriod=5000,command='if O.iter<21000 : yade.qt.center()')
]
#Generate a spiral
@@ -66,7 +66,6 @@
t=O.time)
#yade.qt.Renderer().bound=True
-#O.run(1,True)
plot.plot()
#O.bodies[0].state.angVel=Vector3(0.05,0,0)
Follow ups
