yade-dev team mailing list archive
-
yade-dev team
-
Mailing list archive
-
Message #05195
[Branch ~yade-dev/yade/trunk] Rev 2336: - Some cleaning in ScGeom.
------------------------------------------------------------
revno: 2336
committer: bchareyre <bchareyre@dt-rv020>
branch nick: trunk
timestamp: Thu 2010-07-08 20:14:12 +0200
message:
- Some cleaning in ScGeom.
- Implement precomputed quantities (#define IGCACHE for compiling them).
- Experiment precomputation with ElasticContactLaw (#ifdef IGCACHE), other laws should not be impacted yet whatever IGCACHE.
- Make authors consistent (Galizzi never touched ElasticContactLaw).
modified:
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/ElasticContactLaw.cpp
pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp
--
lp:yade
https://code.launchpad.net/~yade-dev/yade/trunk
Your team Yade developers is subscribed to branch lp:yade.
To unsubscribe from this branch go to https://code.launchpad.net/~yade-dev/yade/trunk/+edit-subscription
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-06-29 14:32:03 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.cpp 2010-07-08 18:14:12 +0000
@@ -6,9 +6,88 @@
#include "ScGeom.hpp"
#include<yade/core/Omega.hpp>
YADE_PLUGIN((ScGeom));
+
// At least one virtual method must be in the .cpp file (!!!)
ScGeom::~ScGeom(){};
+Vector3r& ScGeom::rotate(Vector3r& shearForce){
+ // approximated rotations
+ shearForce -= shearForce.cross(orthonormal_axis);
+ shearForce -= shearForce.cross(twist_axis);
+ //NOTE : make sure it is in the tangent plane? It's never been done before. Is it not adding rounding errors at the same time in fact?...
+ //shearForce -= normal.dot(shearForce)*normal;
+ 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
+ //This line gives nan randomly, when it happens the same repeated lines will give (0,0,0) below
+ orthonormal_axis = prevNormal.cross(normal);
+ // Let's try stupid workaround -> not better
+// Vector3r& lhs= prevNormal;
+// Vector3r& rhs= normal;
+// orthonormal_axis[0]=lhs[1] * rhs[2] - lhs[2] * rhs[1];
+// orthonormal_axis[1]=lhs[2] * rhs[0] - lhs[0] * rhs[2];
+// orthonormal_axis[2]=lhs[0] * rhs[1] - lhs[1] * rhs[0];
+
+ Real angle = dt*0.5*normal.dot(rbp1.angVel + rbp2.angVel);
+ twist_axis = angle*normal;
+ //The previous normal is updated after being used
+// prevNormal=normal;
+ if (isnan(orthonormal_axis[0]) || isnan(orthonormal_axis[1]) ||isnan(isnan(orthonormal_axis[2])))
+ {
+ cerr << "axis "<< twist_axis << " "<<orthonormal_axis << " "<<normal<<" "<<prevNormal<<" "<< rbp1.angVel<<" "<< rbp2.angVel<< endl;
+ cerr<<"nan"<<endl;
+ orthonormal_axis = prevNormal.cross(normal); cerr<<orthonormal_axis<<endl;
+ orthonormal_axis = prevNormal.cross(normal); cerr<<orthonormal_axis<<endl;
+ orthonormal_axis = prevNormal.cross(normal); cerr<<orthonormal_axis<<endl;
+ }
+ //Precompute shear increment
+ Vector3r& x = contactPoint;
+ Vector3r c1x, c2x;
+ if(avoidGranularRatcheting){
+ //FIXME : put the long comment on the wiki and keep only a small abstract and link here.
+ /* B.C. Comment : The following definition of c1x and c2x is to avoid "granular ratcheting". This phenomenon has been introduced to me by S. McNamara in a seminar held in Paris, 2004 (GDR MiDi). The concept is also mentionned in many McNamara, Hermann, Lüding, and co-workers papers (see online : "Discrete element methods for the micro-mechanical investigation of granular ratcheting", R. GarcÃa-Rojo, S. McNamara, H. J. Herrmann, Proceedings ECCOMAS 2004, @ http://www.ica1.uni-stuttgart.de/publications/2004/GMH04/), where it refers to an accumulation of strain in cyclic loadings.
+
+ Unfortunately, published papers tend to focus on the "good" ratcheting, i.e. irreversible deformations due to the intrinsic nature of frictional granular materials, while the talk of McNamara in Paris clearly mentioned a possible "bad" ratcheting, purely linked with the formulation of the contact laws in what he called "molecular dynamics" (i.e. Cundall's model, as opposed to "contact dynamics" from Moreau and Jean).
+
+ Giving a short explanation :
+ The bad ratcheting is best understood considering this small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
+ 1. translation "dx" in the normal direction
+ 2. rotation "a"
+ 3. translation "-dx" (back to initial position)
+ 4. rotation "-a" (back to initial orientation)
+
+ If the branch vector used to define the relative shear in rotation*branch is not constant (typically if it is defined from the vector center->contactPoint like in the "else" below), then the shear displacement at the end of this cycle is not null : rotations a and -a are multiplied by branches of different lengths.
+ It results in a finite contact force at the end of the cycle even though the positions and orientations are unchanged, in total contradiction with the elastic nature of the problem. It could also be seen as an inconsistent energy creation or loss. It is BAD! And given the fact that DEM simulations tend to generate oscillations around equilibrium (damped mass-spring), it can have a significant impact on the evolution of the packings, resulting for instance in slow creep in iterations under constant load.
+ The solution to avoid that is quite simple : use a constant branch vector, like here radius_i*normal.
+ */
+ // 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));
+// cerr <<"relativeVelocity"<<relativeVelocity<<endl;
+ //update relative velocity for interactions across periods
+ //if (scene->cell->isPeriodic) shearDisplacement+=scene->cell->velGrad*scene->cell->Hsize*cellDist; FIXME : we need pointer to scene. For now, everything will be computed outside this function, which doens't make much sense.
+ relativeVelocity+=shiftVel;
+ //keep the shear part only
+ relativeVelocity = relativeVelocity-normal.dot(relativeVelocity)*normal;
+ shearIncrement = relativeVelocity*dt;
+// cerr <<"shearIncrement EOP"<<shearIncrement <<endl;
+}
+
+
+
Vector3r ScGeom::updateShear(const State* rbp1, const State* rbp2, Real dt, bool avoidGranularRatcheting){
Vector3r axis;
@@ -40,11 +119,6 @@
Vector3r c1x, c2x;
if(avoidGranularRatcheting){
- /* The following definition of c1x and c2x is to avoid "granular ratcheting"
- * (see F. ALONSO-MARROQUIN, R. GARCIA-ROJO, H.J. HERRMANN,
- * "Micro-mechanical investigation of granular ratcheting, in Cyclic Behaviour of Soils and Liquefaction Phenomena",
- * ed. T. Triantafyllidis (Balklema, London, 2004), p. 3-10 - and a lot more papers from the same authors) */
-
// FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
c1x = radius1*normal;
c2x = -radius2*normal;
@@ -64,52 +138,22 @@
return shearIncrement;
}
+//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, bool avoidGranularRatcheting){
//Just pass null shift to the periodic version
return rotateAndGetShear(shearForce,prevNormal,rbp1,rbp2,dt,Vector3r::Zero(),avoidGranularRatcheting);
}
+//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){
- Vector3r axis;
- Real angle;
- // approximated rotations
- axis = prevNormal.cross(normal);
- shearForce -= shearForce.cross(axis);
- angle = dt*0.5*normal.dot(rbp1->angVel + rbp2->angVel);
- axis = angle*normal;
- shearForce -= shearForce.cross(axis);
- // exact rotations
- #if 0
- Quaternionr q;
- 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
-
+#ifdef IGCACHE
+ rotate(shearForce);
+ return shearIncrement;
+#else
+ rotate(shearForce,prevNormal,rbp1,rbp2,dt);
Vector3r& x = contactPoint;
Vector3r c1x, c2x;
if(avoidGranularRatcheting){
- //FIXME : put the long comment on the wiki and keep only a small abstract and link here.
- /* B.C. Comment : The following definition of c1x and c2x is to avoid "granular ratcheting". This phenomenon has been introduced to me by S. McNamara in a seminar help in Paris, 2004 (GDR MiDi). The concept is also mentionned in many McNamara, Hermann, Lüding, and co-workers papers (see online : "Discrete element methods for the micro-mechanical investigation of granular ratcheting", R. GarcÃa-Rojo, S. McNamara, H. J. Herrmann, Proceedings ECCOMAS 2004, @ http://www.ica1.uni-stuttgart.de/publications/2004/GMH04/), where it refers to an accumulation of strain in cyclic loadings.
-
- Unfortunately, published papers tend to focus on the "good" ratcheting, i.e. irreversible deformations due to the intrinsic nature of frictional granular materials, while the talk of McNamara in Paris clearly mentioned a possible "bad" ratcheting, purely linked with the formulation of the contact laws in what he called "molecular dynamics" (i.e. Cundall's model, as opposed to "contact dynamics" from Moreau and Jean).
-
- Giving a short explanation :
- The bad ratcheting is best understood considering this small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
- 1. translation "dx" in the normal direction
- 2. rotation "a"
- 3. translation "-dx" (back to initial position)
- 4. rotation "-a" (back to initial orientation)
-
- If the branch vector used to define the relative shear in rotation*branch is not constant (typically if it is defined from the vector center->contactPoint like in the "else" below), then the shear displacement at the end of this cycle is not null : rotations a and -a are multiplied by branches of different lengths.
- It results in a finite contact force at the end of the cycle even though the positions and orientations are unchanged, in total contradiction with the elastic nature of the problem. It could also be seen as an inconsistent energy creation or loss. It is BAD! And given the fact that DEM simulations tend to generate oscillations around equilibrium (damped mass-spring), it can have a significant impact on the evolution of the packings, resulting for instance in slow creep in iterations under constant load.
- The solution to avoid that is quite simple : use a constant branch vector, like here radius_i*normal.
- */
// FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
c1x = radius1*normal;
c2x = -radius2*normal;
@@ -128,9 +172,13 @@
Vector3r shearDisplacement = relativeVelocity*dt;
//shearForce -= ks*shearDisplacement;//this is constitutive behaviour -> moved to Law2 functors
return shearDisplacement;
+#endif //IGCACHE
}
-Vector3r ScGeom::rotate(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt){
+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
@@ -140,40 +188,25 @@
axis = angle*normal;
shearForce -= shearForce.cross(axis);
// exact rotations
- #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
+
+// 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 : put the long comment on the wiki and keep only a small abstract and link here.
- /* B.C. Comment : The following definition of c1x and c2x is to avoid "granular ratcheting". This phenomenon has been introduced to me by S. McNamara in a seminar help in Paris, 2004 (GDR MiDi). The concept is also mentionned in many McNamara, Hermann, Lüding, and co-workers papers (see online : "Discrete element methods for the micro-mechanical investigation of granular ratcheting", R. GarcÃa-Rojo, S. McNamara, H. J. Herrmann, Proceedings ECCOMAS 2004, @ http://www.ica1.uni-stuttgart.de/publications/2004/GMH04/), where it refers to an accumulation of strain in cyclic loadings.
-
- Unfortunately, published papers tend to focus on the "good" ratcheting, i.e. irreversible deformations due to the intrinsic nature of frictional granular materials, while the talk of McNamara in Paris clearly mentioned a possible "bad" ratcheting, purely linked with the formulation of the contact laws in what he called "molecular dynamics" (i.e. Cundall's model, as opposed to "contact dynamics" from Moreau and Jean).
-
- Giving a short explanation :
- The bad ratcheting is best understood considering this small elastic cycle at a contact between two grains : assuming b1 is fixed, impose this displacement to b2 :
- 1. translation "dx" in the normal direction
- 2. rotation "a"
- 3. translation "-dx" (back to initial position)
- 4. rotation "-a" (back to initial orientation)
-
- If the branch vector used to define the relative shear in rotation*branch is not constant (typically if it is defined from the vector center->contactPoint like in the "else" below), then the shear displacement at the end of this cycle is not null : rotations a and -a are multiplied by branches of different lengths.
- It results in a finite contact force at the end of the cycle even though the positions and orientations are unchanged, in total contradiction with the elastic nature of the problem. It could also be seen as an inconsistent energy creation or loss. It is BAD! And given the fact that DEM simulations tend to generate oscillations around equilibrium (damped mass-spring), it can have a significant impact on the evolution of the packings, resulting for instance in slow creep in iterations under constant load.
- The solution to avoid that is quite simple : use a constant branch vector, like here radius_i*normal.
- */
// FIXME: For sphere-facet contact this will give an erroneous value of relative velocity...
c1x = radius1*normal;
c2x = -radius2*normal;
@@ -183,7 +216,6 @@
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
//if (scene->cell->isPeriodic) shearDisplacement+=scene->cell->velGrad*scene->cell->Hsize*cellDist; FIXME : we need pointer to scene. For now, everything will be computed outside this function, which doens't make much sense.
relativeVelocity+=shiftVel;
=== modified file 'pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp'
--- pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-06-21 17:31:31 +0000
+++ pkg/dem/DataClass/InteractionGeometry/ScGeom.hpp 2010-07-08 18:14:12 +0000
@@ -13,13 +13,18 @@
*/
#define SCG_SHEAR
+// #define IGCACHE
class ScGeom: public GenericSpheresContact {
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;
- Real &radius1, &radius2;
-
+ Vector3r shearIncrement;
+
//! 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);
@@ -28,9 +33,17 @@
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 vector as well).
- Vector3r rotate(Vector3r& shearForce, const Vector3r& prevNormal, const State* rbp1, const State* rbp2, Real dt);
+ //! 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);
+
// 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.
@@ -42,7 +55,7 @@
((Vector3r,prevNormal,Vector3r::Zero(),"Normal of the contact in the previous step. |ycomp|"))
,
/* extra initializers */ ((radius1,GenericSpheresContact::refR1)) ((radius2,GenericSpheresContact::refR2)),
- /* ctor */ createIndex();,
+ /* ctor */ createIndex(); shearIncrement=Vector3r::Zero(); twist_axis=orthonormal_axis=Vector3r::Zero();,
/* py */
.def_readwrite("penetrationDepth",&ScGeom::penetrationDepth,"documentation")
);
=== modified file 'pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-06-29 14:32:03 +0000
+++ pkg/dem/Engine/Functor/Ig2_Box_Sphere_ScGeom.cpp 2010-07-08 18:14:12 +0000
@@ -12,7 +12,7 @@
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-common/Box.hpp>
-
+#include<yade/core/Scene.hpp>
#include<yade/lib-base/Math.hpp>
@@ -59,6 +59,7 @@
if (cOnBox_boxLocal[2]<-extents[2]){cOnBox_boxLocal[2]=-extents[2]; inside=false; }
if (cOnBox_boxLocal[2]> extents[2]){cOnBox_boxLocal[2]= extents[2]; inside=false; }
+ shared_ptr<ScGeom> scm;
if (inside){
// sphere center inside box. find largest `cOnBox_boxLocal' value:
// minCBoxDist_index is the coordinate index that minimizes extents[minCBoxDist_index]-abs(cOnBox_boxLocal[minCBoxDist_index] (sphere center closest to box boundary)
@@ -92,19 +93,18 @@
*/
pt1 = se32.position+normal*minCBoxDist;
pt2 = se32.position-normal*s->radius;
-
- shared_ptr<ScGeom> scm;
+ Vector3r normal = pt1-pt2; normal.normalize();
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=pt1-pt2; scm->prevNormal.normalize(); }
+ 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 = pt1-pt2; scm->normal.normalize();
+ scm->normal = normal;
scm->penetrationDepth = (pt1-pt2).norm();
scm->radius1 = s->radius;
scm->radius2 = s->radius;
@@ -139,24 +139,26 @@
pt2=se32.position+cOnBox_sphere*s->radius;
- shared_ptr<ScGeom> scm;
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;}
-
+ 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;
}
+#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);
+#endif
return true;
}
=== modified file 'pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp'
--- pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp 2010-05-04 13:56:05 +0000
+++ pkg/dem/Engine/Functor/Ig2_Facet_Sphere_ScGeom.cpp 2010-07-08 18:14:12 +0000
@@ -10,7 +10,7 @@
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/Sphere.hpp>
#include<yade/pkg-common/Facet.hpp>
-
+#include<yade/core/Scene.hpp>
#include<yade/lib-base/Math.hpp>
CREATE_LOGGER(Ig2_Facet_Sphere_ScGeom);
@@ -95,21 +95,27 @@
if (penetrationDepth>0 || c->isReal())
{
shared_ptr<ScGeom> scm;
+ bool isNew = !c->interactionGeometry;
if (c->interactionGeometry)
scm = YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
else
scm = shared_ptr<ScGeom>(new ScGeom());
normal = facetAxisT*normal; // in global orientation
- scm->contactPoint = se32.position - (sphereRadius-0.5*penetrationDepth)*normal;
+ 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->dt,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-06-29 14:32:03 +0000
+++ pkg/dem/Engine/Functor/Ig2_Sphere_Sphere_ScGeom.cpp 2010-07-08 18:14:12 +0000
@@ -6,7 +6,7 @@
#include"Ig2_Sphere_Sphere_ScGeom.hpp"
#include<yade/pkg-dem/ScGeom.hpp>
#include<yade/pkg-common/Sphere.hpp>
-
+#include<yade/core/Scene.hpp>
#include<yade/lib-base/Math.hpp>
#include<yade/core/Omega.hpp>
@@ -29,19 +29,26 @@
Real penetrationDepthSq=pow(interactionDetectionFactor*(s1->radius+s2->radius),2) - normal.squaredNorm();
if (penetrationDepthSq>0 || c->isReal() || force){
shared_ptr<ScGeom> scm;
- bool isNew=c->interactionGeometry;
- if(c->interactionGeometry) scm=YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
+ bool isNew = !c->interactionGeometry;
+ if(!isNew) scm=YADE_PTR_CAST<ScGeom>(c->interactionGeometry);
else { scm=shared_ptr<ScGeom>(new ScGeom()); c->interactionGeometry=scm; }
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;
+ 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->dt,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;
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-06-16 16:39:17 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.cpp 2010-07-08 18:14:12 +0000
@@ -1,6 +1,5 @@
/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi olivier.galizzi@xxxxxxx *
-* Copyright (C) 2009 by Bruno Chareyre bruno.chareyre@xxxxxxxxxxx *
+* Copyright (C) 2005 by Bruno Chareyre bruno.chareyre@xxxxxxxxxxx *
* *
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
@@ -12,7 +11,6 @@
#include<yade/pkg-dem/DemXDofGeom.hpp>
#include<yade/core/Omega.hpp>
#include<yade/core/Scene.hpp>
-#include<yade/core/Scene.hpp>
YADE_PLUGIN((Law2_ScGeom_FrictPhys_Basic)(Law2_Dem3DofGeom_FrictPhys_Basic)(ElasticContactLaw));
@@ -63,21 +61,23 @@
return;}
State* de1 = Body::byId(id1,ncb)->state.get();
State* de2 = Body::byId(id2,ncb)->state.get();
- Vector3r& shearForce = currentContactPhysics->shearForce;
- Real un=currentContactGeometry->penetrationDepth;
+ Real& un=currentContactGeometry->penetrationDepth;
TRVAR3(currentContactGeometry->penetrationDepth,de1->se3.position,de2->se3.position);
currentContactPhysics->normalForce=currentContactPhysics->kn*std::max(un,(Real) 0)*currentContactGeometry->normal;
- //FIXME : it would make more sense to compute the shift in some Ig2->getShear, using precomputed velGrad*Hsize in Cell and I->cellDist, but it is not possible with current design (Ig doesn't know scene nor I).
+#ifdef IGCACHE
+ Vector3r& shearForce = currentContactGeometry->rotate(currentContactPhysics->shearForce);
+ 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 {
- Vector3r dus = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt,shiftVel,/*avoid ratcheting*/true);
- //Linear elasticity giving "trial" shear force
- shearForce -= currentContactPhysics->ks*dus;
- }
+ } else shearForce -= currentContactPhysics->ks*shearDisp;
// PFC3d SlipModel, is using friction angle. CoulombCriterion
Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
@@ -88,8 +88,6 @@
} else {
//almost the same with 2 additional Vector3r instanciated for energy tracing, duplicated block to make sure there is no cost for the instanciation of the vectors when traceEnergy==false
if(useShear) throw ("energy tracing not defined with useShear==true");
- Vector3r shearDisp = currentContactGeometry->rotateAndGetShear(shearForce,currentContactPhysics->prevNormal,de1,de2,dt,shiftVel,true);
- Vector3r prevForce=shearForce;//store prev force for definition of plastic slip
shearForce -= currentContactPhysics->ks*shearDisp;
Real maxFs = currentContactPhysics->normalForce.squaredNorm()*
std::pow(currentContactPhysics->tangensOfFrictionAngle,2);
@@ -112,6 +110,7 @@
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));
}
+ //FIXME : make sure currentContactPhysics->prevNormal is not used anywhere, remove it from physics and remove this line :
currentContactPhysics->prevNormal = currentContactGeometry->normal;
}
=== modified file 'pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp'
--- pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-06-16 16:39:17 +0000
+++ pkg/dem/Engine/GlobalEngine/ElasticContactLaw.hpp 2010-07-08 18:14:12 +0000
@@ -1,6 +1,5 @@
/*************************************************************************
-* Copyright (C) 2004 by Olivier Galizzi olivier.galizzi@xxxxxxx *
-* Copyright (C) 2009 by Bruno Chareyre bruno.chareyre@xxxxxxxxxxx *
+* Copyright (C) 2005 by Bruno Chareyre bruno.chareyre@xxxxxxxxxxx *
* *
* This program is free software; it is licensed under the terms of the *
* GNU General Public License v2 or later. See file LICENSE for details. *
